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INTRODUCTION 


The  initial  objectives  of  this  study  were  quite  simple, 
namely,  to  analytically  determine  whether  a liquid-filled 
projectile  having  a partial  fill  similar  to  the  XM736,  eight- 
inch  binary  projectile  would  be  an  adequate  ballistic  match 
to  a comparable  solid  projectile,  the  M509.  This  study  is 
an  effort  to  respond  to  a portion  of  a study  request  of  the 
XM736  projectile  initiated  by  the  U.S.  Army  Materiel  Command. 
During  our  investigation  the  objectives  and  techniques  were 
expanded.  To  a great  extent  the  study  assumed  a methodolo- 
gical orientation,  with  the  methods  illustrated  by  simple 
examples . 

In  an  effort  to  keep  the  report  unclassified  some  spe- 
cificity and  fidelity  to  actual  developmental  configurations 
had  to  be  sacrificed.  For  example,  the  effect  of  change  in 
liquid  position  on  inertial  properties  was  examined  via  a 
hypothetical,  simplified,  liquid-filled,  projectile  con- 
figuration. Hopefully,  this  idealization  does  not  preclude 
application  to  real  systems.  All  numerical  examples  given 
here  were  ohoaen  from  the  eight-inch  family  of  ammunition. 
However,  the  methods  used  are  considered  to  be  applicable  to 
other  spin-stabilized,  liquid-filled  projectiles. 

Insofar  as  the  initial  study  objectives  are  concerned, 
one  can  conclude  that  the  XM736  projectile  is  a reasonable 
ballistic  match  to  the  Mf>09.  Difference  in  mean  point  of 
impact  can  be  corrected  by  small  adjustments  in  aiming.  The 
XM73&  is  judged  to  be  ballistically  stable  but  will  likely 
have  a somewhat  larger  ballistic  dispersion  than  either  the 
Mf>09  or  M106  projectiles. 

These  conclusions  certainly  can  be  challenged  on  the 
grounds  that  the  methods  and/or  data  inputs  are  inapplicable 
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to  the  XM736  projectile.  Admittedly,  the  phenomena  occuring 
during  spinup  and  mixing  of  the  chemical  reagents  in  that 
system  are  not  treated  adequately  in  this  study.  In  fact, 
limitations  of  the  study  are  recognized  and  stated  explicitly 
in  the  caveats  below. 


Although  lacking  physical  rigor  at  several  points, 
this  report  is  offered  in  the  hope  that  some  of  the  methods 
and  insights  may  be  of  interest  and  that  some  of  the  physical 
shortcomings  may  provoke  further  inquiry.  The  author  suggests 
that  a careful,  finite-element  approach  to  the  analysis  of 
mixing  and  spinup,  considering  pertinent  physical  phenomena, 
may  be  a useful  avenue  of  approach. 


Caveats 

The  following,  numbered  caveats  are  provided  to  iden- 
tify restrictions  in  the  scope  and  depth  of  the  study: 

1.  Treatment  of  phenomena  occuring  during  mixing  of 
chemical  reagents  (for  binary  systems)  is  highly  idealized. 

a.  Constant  free  volume  within  the  shell  cavity  was 
assumed  although  it  is  known  that  prior  to  mixing 
the  liquid  ingredients  in  the  XM736  occupy  a smaller 
volume  than  that  assumed,  and  that  after  mixing  and 
reacting  chemically,  essentially  the  entire  cavity 
is  occupied  by  reaction  products. 

b.  The  value  of  kinematic  viscosity  was  treated  parame- 
trically at  two  (constant)  levels  --  10  and  1 centistoke. 
It  is  estimated  from  laboratory  tests  that  the  viscosity 
of  the  cold  reactants  is  greater  than  ten  centistokes 
whereas  the  hot  reaction  products  have  a viscosity  of 
less  than  one  centistoke  • 

2.  Only  certain  kinds  of  liquid  resonances  were  treated. 

It  is  known  that  pressure  waves  propagating  in  a liquid 
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confined  within  a cavity  have  certain  characteristic-  or 
eigen-frequencies.  Although  not  treated  in  this  report,  a 
discussion  of  these  resonant  frequencies  is  found  in  AMCP- 
706-165,  Liquid-Filled  Projectile  Design.  The  discussion  of 
liquid  resonances  in  this  report  is  limited  to  the  whole- 
body  or  sloshing  motion  of  the  liquid.  Obviously,  this 
discussion  applies  only  to  partially-filled  shells. 

3.  An  academic  or  idealized  mathematical  treatment  of  liquid 
spinup  is  contained  herein.  This  development  strictly  applies 
only  when  liquid  flow  is  laminar  and,  then,  only  for  cavities 
which  are  quite  long  relative  to  their  diameter.  The  dif- 
fusion of  angular  momentum  within  the  liquid  due  to  turbulence, 
generated  either  by  a chemical  reaction  or  by  vorticity  cells 
induced  by  the  boundary  conditions  at  the  ends  of  the  liquid 
cylinder,  were  not  considered  in  this  study.  For  most 
liquid-filled  projectiles,  the  rapid  achievement  of  homo- 
genous angular  velocity  is  due  mainly  to  the  turbulent  con- 
dition of  the  liquid  during  spinup  which,  effectively,  creates 
a large  apparent  viscosity. 


Next  page  is  blank. 
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CHAPTER  1 

CHANGE  IN  INERTIAL  CHARACTERISTICS  OP 
A LIQUID-PILLED  PROJECTILE  DURING  PLIGHT 

Inertial  Characteristics  In  Two  Limiting  Configurations 

Consider  the  following  idealized  model  of  a liquid- 
filled  projectile  during  launch.  Because  of  the  large  axial 
acceleration  and  small  average  angular  velocity  of  the  liquid, 
the  liquid  will  be  forced  to  the  rear  of  the  shell  cavity  so 
that  its  free  surface  will  approximate  a flat  circular  dia- 
phragm. See  Example  1 below.  This  liquid  configuration.  A, 
is  shown  below  in  Figure  1.1. 


The  solid  shell  body  is  approximated  as  a cylindrical 
sleeve  32  in.  long,  with  a 1 in.  wall  thickness.  Massless 
boundaries  confine  the  liquid  to  a cylindrical  cavity  22  in. 
long,  having  a 6 in.  diameter.  The  liquid  fills  only  90 % of 
the  available  volume.  During  launch  the  height  of  the  column 
of  liquid  is  19.6  in.  in  the  axial  or  X-direction.  The  average 
density  of  the  solid  shell  body  is  assumed  to  be  0.2558  lb/in^, 
i.e.,  specific  gravity  7*09.  The  density  of  the  liquid  is 
assumed  to  be  0.0361  lb/in  , i.e.,  specific  gravity  1.  See 
ref.  [1]*  for  standard  formulas  and  material  values. 

# Square-bracketed  numbers  refer  to  the  bibliographic  citations. 
[1]  Eshbach,  C.  Handbook  of  Engineering  Fundamentals,  c.  1952. 
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Other  pertinent  parameter  values  are  shown  in  Table  1.1. 

TABLE  1.1.  INERTIAL  CHARACTERISTICS  OF 
A HYPOTHETICAL  LIQUID-FILLED  PROJECTILE  IN  CONFIGURATION  A 


Characteristic 

Value 

Dimension 

volume  of  solid 

703.7 

in3 

weight  of  solid 

180 

lb 

volume  of  liquid 

554-2 

in3 

weight  of  liquid 

20 

lb 

total  weight  of  projectile 

200 

lb 

center  of  gravity  of  liquid,  x^ 

9.8 

in 

center  of  gravity  of  solid,  x^ 

16 

in 

center  of  gravity  of  projectile 

15.38 

in 

axial  moment  of  inertia  of 
projectile 

0.50506 

slug 

ft2 

axial  moment  of  liquid 

0.019426 

slug 

ft2 

axial  moment  of  solid 

0.48564 

slug 

ft2 

pitch  moment  of  inertia  of  pro- 
jectile about  its  eg 

3.8554 

slug 

ft2 

pitch  moment  of  liquid  about  x^ 

0.14791 

slug 

ft2 

pitch  moment  of  solid  about  x^ 

3.55813 

slug 

ft2 
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After  exit  from  the  cannon,  the  projectile  will  ex- 
perience negligible  axial  acceleration,  relative  to  its 
centripetal  acceleration.  Therefore,  the  free  surface  of 
the  liquid  will  be  a cylinder  centered  on  the  spin  axis. 
This  configuration,  B,  is  shown  in  Figure  1.2. 


In  Configuration  B,  the  volumes  and  weights  of  solid 
and  liquid  are,  of  course,  the  same  as  in  Conf iguration  A. 
However,  other  inertial  properties  have  changed  from  those 
shown  in  Table  1.1  to  the  values  shown  below: 

center  of  gravity  (eg)  of  liquid  s 11.0  in 

center  of  gravity  (eg)  of  projectile  ■ 15-5  in 

Thus,  the  projectile  eg  has  moved  forward  by  0.12  in.  In  the 
M10&  projectile,  which  has  the  same  mass  and  caliber  as  this 
hypothetical  projectile,  the  round-to-round  standard  devia- 
tion of  eg  position  about  the  mean  is  0.06lj.  in.,  due  to 
manuf acturing  tolerances  [ 2} . Consequently,  one  may  conclude 
that  eg  shifts  due  to  liquid  configuration  changes  in  this 

[2]  Final  Report  of  the  HY-BRA  Weapons  System,  Feb.  1970. 
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hypothetical  projectile  are  less  than  two  times  the  standard 
deviation  associated  with  round-to-round  eg  shifts  in  com- 
parable solid  projectiles.  Also 

axial  moment  of  inertia  of  liquid 

(Configuration  B)  ■ 0.02151*.  slug  ft2 

axial  moment  of  inertia  of  projectile 

(Configuration  B)  * 0.50718  slug  ft2  . 

Further,  comparisons  of  the  tabular  data  show  that  the 

axial  moment  of  inertia  for  Configuration  B is  0.00212  slug 

p 

ft  greater  than  for  Conf iguration  A.  To  put  this  difference 
in  perspective,  consider  that  the  estimated  round-to-round 
standard  deviation  in  spin  inertia  for  the  M106  projectile, 

O 

due  to  manufacturing  variability,  is  0.00252  slug  ft  [2]. 
Thus,  the  change  in  spin  inertia  resulting  from  liquid  con- 
figuration changes  is  less  than  one  standard  deviation  of 
that  of  the  comparable  solid  projectile.  Using  the  deflec- 
tion sensitivity  to  spin  inertia  given  in  Reference  [2]  for 
the  M106  projectile,  one  expects  a change  in  deflection  of 
0.26  mils  due  to  change  in  inertial  properties  associated 
with  the  change  in  configuration. 

Further, 

pitch  moment  of  inertia  of  liquid 

(Configuration  B)  about  x^  - O.I8I4.88  slug  ft2 

pitch  moment  of  inertia  of  projectile 

p 

(Configuration  B)  about  its  eg  = 3»81j.01  slug  ft 


[2]  Ibid. 


Note  that  even  though  the  pitch  inertia  of  the  liquid 
d>out  its  eg  is  greater  in  Configuration  B than  in  Configu- 
ration A,  the  pitch  moment  of  inertia  of  the  whole  projectile 
is  smaller  for  Configuration  B.  This  is  due  to  the  shift 
in  eg  of  the  liquid  which  brings  it  closer  to  the  eg  of 

the  whole  projectile.  The  difference  in  pitch  inertia 

2 

between  the  oonf igurations  is  0.0152  slug  ft  . 

The  estimated  [2]  round-to-round  standard  deviation  of 
pitch  inertia  in  the  M106  projectile  due  to  typical  manu- 

O 

facturing  tolerances  is  0.0252  slug  ft  . Thus,  liquid 
configuration  change  in  the  hypothetical  projectile  pro- 
duces a change  in  pitch  inertia  only  0.6  of  one  standard 
deviation  in  pitch  inertia  of  the  M106. 

One  must  conclude  that  the  changes  in  inertial  charac- 
teristics accon^anying  liquid  conf igurational  change  in  the 
hypothetical  liquid-filled  projectile  would  not  be  suffi- 
cient to  produce  significant  changes  in  exterior  ballistics, 
given  only  adequate  ballistic  stability  of  the  liquid-filled 
projectile. 


[2]  Ibid 


Estimate  of  Shape  of  the  Free  Surface  of  the  Liquid 
In  a Liquid-Filled  Projectile  During  Acceleration 


Consider  Figure  1.3  below  in  which  a cross  section  of 
the  liquid-filled  cavity  of  a projectile  is  depicted.  Ac- 
celeration of  the  projectile  is  assumed  in  the  positive 
X-direction. 


Y 

A 


m 


x V,  free 

surface 


Figure  1.3*  Liquid  Free  Surface  During  Acceleration 


The  free  surface  is  defined  by  the  function  y(x). 

The  differential  element  at  (x,y),  dv,  experiences  forces 
in  the  Y-  and  X-directions  due  to  spin  and  setback,  re- 
spectively. Neglecting  gravity,  which  is  small  compared 
to  the  accelerations  of  interest,  and  assuming  solid-body 
rotation,  the  forces  on  dv  may  be  written  as 

Fx  = " ax  p dv 

Py  a P dv  y to2  (1.1) 

with 

dv  s 2 rr  y dy  dx 
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The  angular  velocity, 
velocity  of  projectile*  v , 


o),  is  given  in  terms  of  the 
caliber,  D,  and  twist,  T,  as 


0) 


2 It  V 

E 


D T 


(1.2) 


The  vector  sum  of  these  forces  on  dv  must  be  normal  to 
the  free  liquid  surface  at  equilibrium, since  the  liquid 
will  not  support  shear  stresses.  Therefore,  the  slope  of 
the  free  surface  at  (x,y)  is  given  by 


• - (tan  3)'1 


-F 

~ T 


X , -1 

= k y 


with 


k = a*  “ 


-2 


(1.3) 


Equation  (1.3)  may  be  integrated  to  yield 


y2  = 2 k (X  - xc) 


or 


y - [2  k (x  - xo)3 


X, 

'2 


(14) 


and 


x * xo  + ^ ■ 


(1.5) 


where  the  value  of  xq  is  determined  by  the  volume  of  the 
cavity  and  volume  of  liquid.  Since  the  volume,  V,  of  liquid 
is  assumed  constant. 


p y i 

V " 71  (ym  ~ yi)  xm  = J 2 % y x(y)  dy’ 


(1.6) 


0 - ^ y ,0<x  ^ x 

Jjw  9 n m » 


m 
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with 


(1.7) 


7l  = minty,,,  y(x„)  = ( 2 k (x^  - mD)  J ) . 

Then,  from  (1.5) 


1 3 


¥ * <4  - T-  = / ( A + *o  y)  ** 


JL  yrn  % 

2 TT  ' 2 


2 ll 

7i  y} 


^-^r-rr-0 


For  yx  < ym. 


4 ■ ¥ <*.  - v> 


(1.8) 


with 


V s it  x y* 
c m Jm 


And  for  - y , 


xo  s 


w y, 


y3 

Jm 

IP?  * 


m 


(1.9) 


Example  1. 

Take  the  following  conditions  at  peak  projectile  ac- 
celeration: 


liquid  volume,  V - 553*8  in^  (0.3205  ft^) 

cavity  volume,  V ■ 622  in^  (0.36  ft^) 

c 

projectile  velocity,  v^  - 200  ft/s 

gun  tube  twist,  T ■ 20  oal/rev 

caliber,  D = 8 in  (2/3  ft) 

projectile  acceleration,  a^  - i|..504  10^  ft/s 

x^  s 22  in  (1.833  ft) 

ym  a 3 in  (0.25  ft)  . 
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2 n v 

s = 94- - 2L|_8  rad/sec  (15  hz)  , 

D T 

a 

= — = 50.70  ft  , 
to2 

y2  = 1014  (x  - xQ) 

Assuming 

*1  = Jm  ■ °'2S  ft- 
xty^)  - xc  - 7 10'3  in 

and 

xq  = 19. in. 


Then , 

co 


k 
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Y 


A 


C>  X 


Pigura  1.4.  Liquid  Free  Surface  During  Acceleration 


Figure  I.I4.  represents  the  position  of  the  liquid  sur- 
face where  the  unoccupied  volume  of  the  cavity  is  large. 

In  this  oase  the  free  (unoccupied)  volume  is 


^f  = ^cavity  " 1 liq 


(1.10) 


From  Figure  1*4» 

r x, 

Vf  = 71  ^xm  “ xl^  ym  + J 71  y2  dx,  0 < x < xm  (1.11) 

xo 


Using  (lv4,  1.11), 

Vt  * ,(lm  - xl>  4 + ” k(xl  ‘ xo)a 

But,  from  (l.f>). 


(1.12) 


(1.13) 
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Then.  (1.12,  1.13)  yield  an  expression  for  x^. 


X1  = xm 


(1.11*) 


Given  the  cavity  dimensions  and  liquid  volume,  Equation 

(l.llj.)  can  be  solved  for  x^  and  with  this  value  xq  can  be 
obtained  from  (1.13). 

Example  2 

Suppose  that  conditions  at  the  muzzle  of  a gun  are  as 
follows 

projectile  acceleration, a - 8.086  10^  ft/s 

X 

projectile  velocity,  v^  = 2000  ft/s 
twist, T - 20  cal/rev 
caliber,  D * 2/3  ft  (8  in) 
cavity  free  volume,  = 0.03935  ft^ 
cavity  length,  x^  - 1.8333  ft 
cavity  radius,  y^  = 0.25  ft 
Using  (1.2), 

w = 914.2.4.8  rad  sec-'*'  ( 150  hz) 

Prom  (1.3), 

k = 0.09103  ft  , 
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and 


y2  = 0.18206  (x  - xq)  (ft)2 
Then  ( 1 . ll|.)  yields 

x±  = 1.80ij.57  ft  = 21.655  in  , 
and  from  (1.13), 

XQ  = 1 .ij.6128  ft  = 17.535  in  , 
so  that 

x,  - x - k.12  in. 

1 o ^ 

The  radial  acceleration  at  y in  this  example  is  a . 

Jm  y 

ay  (Jm)  = ym  = • 25  ( 9ij.2 . ij.8 ) 2 

ay(ym}  = 22*2  1Q,k  t/sZ 
= 0.69  10^  g . 
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CHAPTER  II 

ANGULAR  ACCELERATION  OP  THE  LIQUID  IN  A 
LIQUID-PILLED  PROJECTILE  DURING  PLIGHT 
Liquid  "Splnup”  In  Configuration  A 

We  treat  the  position  of  the  liquid  within  the  shell 
as  in  Configuration  A (Figure  1.1).  Let  the  tangential  or 
circumferential  velocity  of  the  liquid  at  radius  r and  time 
t be  denoted  by  v.  We  neglect  axial  velocity  components. 

Consider  the  annular  volume  element  shown  in  Figure  2.1. 


Figure  2.1.  Volume  Element  Under  Angular  Acceleration 

At  position  (r,x)  in  the  liquid,  the  shear  stress 
acting  on  the  surface  dA_  (below  r)  is 


dA  = 2 n r dx 

At  a radial  position  r + dr  the  shear  stress  acting  on 
the  surface  dA  is 


- or(r,x) 


where 


or(r+dr,x)  , 


-13- 


where 


dA  = 2 ti  (r+dr)  dx 
+ 

Similarly  at  the  left  of  the  element  the  shear  stress 
acting  on  the  differential  area  dAQ  is 

- ax(r,x) 


where 

dA  = 2 tc  r dr 
o 

And  at  axial  position  x + dx,  the  shear  stress  acting 

on  dA  is 

0 

ox(r,x+dx) 

The  axial  moment  of  inertia  of  the  differential  element 
i s 

1 = 2 it  p r^  dx  dr  . (2.1) 

The  sum  of  all  torques  acting  on  the  element  is 
Mj  = [ox(r,x+dx)  - ox(r,x)]  dAQ  r 

+ [or(r+dr,x)  dA+  - or(r,x)  dA_]  r . (2.2) 

Defining  the  angular  velocity  of  the  liquid  at  (r,x)  at 
time  t as  o)(r,x,t),  Newton's  law  gives 

ro  = VLj/1  , (2.3) 
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where  the  functional  dependence  on  r,x  and  time  t has  been 
suppressed  notationally . 

Then 


co  = (Pr)"1  [a  (r,x+dx)  - o (r,x)J  / dx 

A A 

+ (pr2)"1  [ or( r+dr ,x)  ( r+dr ) - or(r,x)  r]  / dr 
Or  in  the  limit  as  dx  and  dr  approach  zero, 

-1  -1  a(r0J 

o)  = (pr)  1 — £■  + (pr2)  1 — — — . ( 2 . Ip) 


But  the  dynamic  viscosity  n is  defined  such  that 


o 

X 


T) 


r 


3m 

3x 


a 


r 


T) 


(2.5) 


Then  with  T]  taken  independent  of  x,r 

co  = (pr)"1  T)  r + (pr2)"1  T)  — (r2  |&)  . 

w = V^  + i>r~2  [2r  g + r2  f-y] 


where  U is  the  kinematic  viscosity. 
V - B p"1 


(2.6a) 


co 


= V 


3 3 co 

3xS 


+ 


2 r 


-1 


Bco 

Br 


B2co*| 

bP^J 


(2.6b) 
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If  the  axial  variability  in  go  can  be  neglected,  the 
simplified  result  in  one  spacial  dimension  is 


The  tangential  velocity  v in  the  liquid  at  radius  r is 


where  the  subscripts  indicate  partial  differentiation  with 
respect  to  that  variable.  Note  that  for  stationary  condi- 
tions, i.e.,  zero  w^,  wf  is  a constant  c and  w(r)  = c r. 

Notation  can  be  simplified  by  defining  a dimensionless 
radius  e and  dimensionless  time  t . 

Let 


(2.7) 


v = r co 


(2.8) 


with  v a function  of  r and  time  t 


v 


w(r, t) 


(2.9) 


With  (2.7,  2.8) 


or 


w 


t 


rr 


0 , 


(2.10) 


€ = r/rx 


and 


t = t V / r 


1 


(2.11) 
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with  r^,  V constants. 

Then  Equation  (2.10)  becomes 


3 u 

3 T = 


u 


9 C 


0 < % < 1 
T > 0 


where 

u(^,t)  = w(r,t) 

The  boundary  conditions  for  (2.12)  are 
u(£,0)  =0  , 0 $ ^ ^ 1 
u(0,t)  =0  , t > 0 
u(  1,t)  = <p(t) 

with 

<P(t)  = f(t)  , t > 0 

f(t)  = fx(t)  + f2(t) 

f^(t)  = a t , t^O 

f2(t)  =0  , t s t± 

f^(t)  — — a(t  — t-^)  > t > t^  • 

Thus 

'P(t)  = + ^2 


(2.12) 


(2.13) 


(2.  111.) 
(2.15) 


(2.16) 


(2.1?) 
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•P^(t)  = *<  T , T ^ 0 

<Pp(T)  =0  , T ^ T1 

'■P2(t)  r - -<(t  - Tx)  , T > T-j^  , 

with  * 

T1  * ^1  ^ / ri 

Taking  the  Laplace  transform  of  (2.12) 
s u”(£,  s)  - u(p,  0)  = u“  (*f  s) 
or  from  ( 2 . 1L|. ) 

s vT(r,  s)  = uJT  (5,  s)  . 

Also,from  (2.15,  2.16), 
u *”  ( 0 , s ) = 0 , 

and 

( 1 , s ) = (f  " ( s ) , 


or 


q>  ” ( s ) - x s 


-2 


- -< 


-2 

s e 


-Tls 


(2.18) 


(2.19) 


( 2 . 20  ) 


(2.21) 


(2.22) 


Since  u“(^,  s)  is  not  a function  of  time,  the  partial 
derivatives  in  Equation  (2.20)  are  actually  total  deriva- 
tives  of  u with  respect  to 
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Thus  , 


da  u'~ 
d r2 


s u” 


= 0 


(2.23) 


whose  solution  with  initial  conditions  given  by  (2.21)  is 

sinh  fs  t. 


u ” ( £ » s ) = <p  " ( s ) 


sinh  /s' 


(2.21+) 


Using  a series  expansion  of 
sinh  /a  p / sinh  'Js1 


--  see  p.  139,  Churchill  [3]  --  and  applying  the  convolution 
theorem 


n=o 


<p(x  - (2  n + 1 - £)*  / 0+  ^))  e-^  d\ 


(2  n + 1 + p)z  / (1+  \*)) 


e 


-\2 


d\ } 


t 


with 

^(n,5) 


2 n + 1 - g 

2 \fx 


^2(n>?) 


2 n + 1 + % 

2 -fT 


(2.25) 


[33  Churchill,  R.V.  Operational  Mathematics,  2nd  Ed., 
McGraw-Hill,  New  York,  c.  195&* 
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For  the  special  case  of  <p(t)  a constant  cp 


o 


2 n + 1 


2/F 


cr 


(2.26) 


Substitution  of  (2.22)  into  (2.25)  appears  to  be  too 
complex  to  pursue.  A numerical  approach  starting  with 
Equation  (2.10)  and  boundary  conditions  given  by  (2.11;, 

2.15»  2.16)  was  judged  more  profitable. 

For  low  viscosity  liquids  and  rapid  rise  of  driving 
angular  velocity  to  a nearly  constant  level*  the  result  in 
(2.26)  may  be  a reasonable  description.  This  expression 
was  evaluated  for  unity  <pQ  and  plotted  in  Figure  2.2. 

Generally  it  was  possible  to  truncate  the  series  at  four 
terms  or  less  to  achieve  0.1$  accuracy. 

The  tangential  velocity  of  the  liquid  in  Configuration  A 
is  also  shown  as  a function  of  t for  several  radial  positions 

in  Figure  2.3. 

Because  angular  acceleration  of  the  liquid  varies  with 
position,  it  is  useful  to  define  an  effective  angular  velocity 
as  that  value  which  would  give  the  liquid  angular  momentum  if 
distributed  uniformly  thruout  the  liquid. 

Thus,  for  the  nondimens ional  case, 


o 


(2.27) 


A plot  of  versus  t is  shown  in  Figure  2.i;. 
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Figure  2.2.  Tangential  Velocity  of  Liquid  Versus  Radial  Position  at  Several  Values  of  Time 
(Liquid  Configuration  A) 


u(£,t) 

2 


Figure  2.3.  Tangential  Velocity  of  Liquid  Versus  Time  for  Three 
Radial  Positions 


22 


Figure  2.4. 
Versus  Time 


Effective  Angular  Velocity  of  Liquid  in  Configuration  A 
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Figure  2.3  shows  that  a nondimens ional  time  % of  0.1  is 

required  to  achieve  50$  or  more  of  the  asymptotic  velocity 
for  E,  > 0.5*  To  apply  these  results  to  a specific  example, 
suppose  that 

V = 10  centistokes  (0.1  stoke) 
r^  = 7*62  cm  (3  in) 

Then  if  "spinup"  is  arbitrarily  defined  as  the  time  at  which 
the  angular  kinetic  energy  reaches  70$  of  its  asymptotic 
value”,  a time  corresponding  to  t = 0.125  is  required. 

This  time 

t = t r^  LT1 

= 0.125  (7.62)2  10 

t = 72.6  sec 

Since  this  time  exceeds  the  maximum  time  of  flight,  one 
would  not  expect  spinup  to  occur  in  a liquid-filled  pro- 
jectile of  this  type  having  the  above  viscosity  and  re- 

•5HC- 

maining  in  Configuration  A. 

Further  during  the  inbore  period,  taken  as  lip  milli- 
seconds, essentially  no  rotation  of  the  liquid  occurs  since 
this  time  corresponds  to  t of  2.1;  10  J for  the  above  example. 

In  order  for  the  liquid  to  remain  in  a stable  Configura- 
tion B after  launch,  the  minimum  radial  acceleration  must 

-x-  Approximately  75$  of  asymptotic  angular  momentum. 

•shs-As  indicated  in  the  Introduction,  the  conditions  for  this  be- 
havior do  not  obtain  in  the  XM73&  projectile,  where  experimen- 
tal data  indicate  that  spinup  is  complete  within  one  to  five 
seconds . 


be  at  least  one  g.  Thus,  the  limiting  angular  velocity  of 
liquid  for  Configuration  B is  19-74  rad/sec;  and  for  a case 
spin  of  942.5  rad/sec,  implies  a limiting  value  of  w.  of 
approximately  0.02.  Frost  Figure  2.4  one  notes  that  this 
value  of  is  obtained  at 

t = 0.012  , 

which,  for  this  example,  corresponds  to  a time 

t = t r£  If  "*■ 

t = 0.012  (7-62)3  10 

t = 7 sec  . 

Beyond  the  time  at  which  Configuration  B is  attained, 
the  spin  dynamics  of  the  liquid  must  be  obtained  from  the 
solution  of  the  following  boundary-value  problem. 

Angular  Acceleration  of  the  Liquid  in  Configuration  B 

As  previously  developed,  the  governing  differential 
equation  is 

"t  =v  [“xx  + 2 r"1  “r  * “rr]  • 

Letting 

u(£,\,t)  = co(r,x,t)  , 


with 
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un 

o 

II 

o 

A.  = x/r1  , 

A.,  = x /r,  , and 

1 ml 

t = t V/v\  , (2.28) 

\ = u\\  + 2 S"1  ur  + » -o  < ^ < 1 (2.29) 

0 < X < \1 

T > 0 

The  initial  and  boundary  conditions  in  this  case  are 


u(§,A.»0)  = 

♦(€»*)  > %Q  < % < 1 

u§(50A,t) 

= o 0 < X < \ 

u(1,A.,t)  = 

<P(t)  0 ^ A.  ^ A.^ 

o 

1! 

<P(t) 

u(^,\-l»t)  = 

= 'P(t)  , T > o (2.30) 

Note  t equal  to  zero  corresponds  to  the  time  at  which  Con- 
figuration B is  realized. 

We  particularize  the  functions  \j/(^)  and  <p(t)  as  followsr 

At  a free  surface  the  shear  stress  in  the  liquid  must 
vanish. 
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(2.31) 


<p(t)  = 1 (Case  1)  , or. 


cp(tt)  = (1.2+  ct) 


-0.2808 


(Case  2)  , 


with 


to.  . and  c constants  . 
lim 


(2.32) 


The  limiting  or  minimum  angular  velocity  for  Configuration 
B is  The  form  of  cp(t)  for  Case  2 is  motivated  by  the 

despin  characteristics  of  a typical  eight-inch  projectile. 

A generalization  of  this  formula  is  developed  in  the  Appendix. 

In  the  case  considered 


r^  = 7*82  cm 
v - 0.1  and  0.01  stoke 

It  is  convenient  to  obtain  a numerical  solution  to  this 
problem  thru  spacial  discretization. 

The  interval  in  £ (^G>1)  is  divided  into  n-1  equal  seg- 
ments of  length 


c = 0.0288  r'£  y'1 


(2.33) 


A?  = (1-S0)  / (n-1) 


(2.34) 
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Similarly,  the  interval  in  K (0,^)  is  divided  into  m-1  equal 
segments  of  length 

= \^  / (m-1)  . (2.35) 

We  define  the  value  of  u at  the  nodal  points  of  this 
segmented  space  as  follows 

uji  = = u^j,Xi»t)  » 1 c ^ n (2.36) 

1 s i m , 

with 

- ?0  + (J-l)  A? 

= (i-1)  A X • 

Division  of  the  space  is  illustrated  in  Figure  2.5. 


u 


n,m 


u 


n,i 


u 


n,  1 


Figure  2.5*  Nodal  Points  for  Finite-Difference  Method 
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The  second  central-difference  approximation  is  used 
to  approximate  u^  and  u^  as  follows: 

- <A?)'2  <Uj+l,i  • 2 uj,i  + V,!1  ’ 


-2 


"U(5j’ki>  - (aX)  (uj,i+l  - 2 uj,i  + 


1 < j < n 
1 < i < m 


(2.37) 


A central-difference  approximation  is  also  used  for  u 


u , , * — u * i . 

u(?  , 
C J 2 A ? 


(2.38) 


1 < j < n 
1 < i < m 


With  these  approximations,  Equation  (2.29)  is  approxi- 
mated as 


dr  uj,i  " <A?)  2 (1*A’/=J)  - ? uj,i 


+ d+ic/sj)  u3+i,ij  + 


(AM  tuj,i-l  " 2 uj,i  + u j , i+ 1 J ’ 


1 < j < n 

1 < i < m . (2.39) 


-29- 


At  the  boundaries  of  the  discrete  space,  from  (2.31,  2.32) 


a?  uj,i  = •'t 


cp  = - 0.2808  c (1.2  + c t) 


-1.2808 


-s—  u . =9  , 1 ^ i 5 n , 

dT  j,m  t ’ J * 


and, 


T~  U 
dT  n,  1 


\ ' 


1 « i =£  m 


(2.40) 


For  Case  1,  c is  set  to  zero. 


- (A?> 


-2  ^ + U 


(u2,i  - ul,i> 


+ (lM"2  1#1,M  '2V+"l,W1  ’ 


1 < i < m 


(2. 4D 


The  first  term  on  the  r.h.s.  of  (2.41)  derives  from  the 
fact  that  the  net  torque  due  to  radial  gradients  acting  on 
the  mass  element  within  the  annular  segment  ( r^  $ r < r~) 
produces  an  angular  acceleration 


W1 


(r2  oa  + r,  a, ) 
2 a r 


9 


and,  for 


(^2  - W-, ) 

= t]  r8  , 

A r 
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(fOe  — co^) 


* 


A r 

[ 1 + (rc/r^)3] 


= -v  (Ar)  2 


2 


((■'a  - w^) 


This  expression  is  equivalent  to  the  first  term  of  the  r.h.s. 
of  ( 2. [j.1 ) . 

The  initial  conditions  for  Cases  1 and  2 (2.31,  2.33)  are 


Note  that  the  initial  projectile  spin  is  normalized  to  unity. 
The  actual  liquid  spin  may  be  obtained  from  the  relation 


figuration  B is  achieved,  the  normalized  spin  is  taken  as 
0 . 9I4J+.  This  time  is  approximately  7 seconds  after  launch 
in  the  examples  treated  here. 

Because  of  the  relatively  weak  longitudinal 

gradients  in  the  angular  velocity  relative  to  the  radial 
gradients,  the  longitudinal  (\)  dependence  was  neglected  in 
the  numerical  results  presented  here. 

With  no  dependence  on  the  index  i.  Equation  (2.39) 
can  be  simplified  as 

* Check  mins  considering  ^-dependence  show  essentially  equi- 
valent results. 


un,i  = 1»  1 « i - m 


u.  . = = 0.02  , 1 < j < n 

j,i  lim  J 


1 < i < m 


(2.^2) 


co(t)  = ooQ  u(e;  ,\  ,t) 


(2.43) 


where  coQ  is  the  projectile  spin  at  launch.  At  the  time  Con- 
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(2.1+4) 


= (A?)-2  [(l-A5/€j)  uj_1  - 2 ^ 
+ (l+A?/?j)  uj+1)  » 1 < j < n 


And,  from  (2.1+0,  2.1+1) 


d u 


n 


dx 


= 0 ( Case  1) 


d u 


= U?)"2  (us  - ux) 


(2.1+5) 


Using  values  of 


n = 11 


€o  = °*33  ’ 

A£  = 0.067  , 


and  a finite  time  step  At  = 10-^,  (2.1+1+,  2.1+5)  were  inte- 
grated using  a fourth-order  Runge-Kutta  procedure.  The 
numerical  results  are  illustrated  in  Figure  2.6. 

The  equivalent  angular  velocity  of  the  liquid  can  also 
be  obtained  for  Configuration  B.  This  is  obtained  numerically 
as  follows : 


I 


n-j. 


J-l 


(?j  + 5M)(uj  + 


WB  = 


(2.1+6) 


j=i 


n-l 


+ "j+1 


This  result  is  plotted  in  Figure  2.7* 


-32- 


T 


Figure  2.6.  Angular  Velocity  of  Liquid  Versus  Time  for  Several 
Radial  Positions  (Liquid  Configuration  B) 
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T 


Figure  2.7*  Effective  Angular  Velocity  of  Liquid  in 
Configuration  B Versus  Time 
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Figure  2.8.  Spin  Versus  Time  for  a Modified  M509  Projectile  at  Zone  7 from 
the  XM201  Cannon  (Inertial  Characteristics  of  the  XM736  Projectile) 


Based  upon  the  run  down  of  spin  over  time  for  an  M$09 
projectile  type,  a second  numerical  example  was  examined 
(Case  2).  The  spin  versus  time  for  a modified  M509  pro- 
jectile is  shown  in  Figure  2.8.  The  boundary  conditions 
at  for  Case  2 — Equation  (2.i|.0)  — closely  match  the 
spin  characteristics  of  the  M509.  Consequently,  this 
example  is  a reasonable  simulation  of  the  spin  dynamics  of 
the  liquid  in  the  XM736  projectile.  As  indicated  in  Equation 
(2.33),  two  values  of  liquid  viscosity  were  used  for  Case  2 — 
1 and  10  centistokes.  Because  of  the  strong  dependence  of 
viscosity  upon  temperature,  this  viscosity  range  is  necessary 
to  encompass  the  expected  range  of  liquid  temperature . In 
this  example  a launch  spin  a>  of  923  rad/sec  is  assumed. 
Numerical  results  are  shown  in  Figure  2.9. 

At  this  point  it  is  necessary  to  repeat  the  caveat  that 
the  above  analysis  of  the  dynamics  of  liquid  spinup  applies, 
strictly,  only  to  liquid-filled  projectiles  having  long, 
narrow  cavities  and  quite  viscous  fills  so  that  laminar  flow 
obtains.  In  the  XM73&,  these  conditions  do  not  hold  (even  at 
10  centistokes).  However,  due  to  the  brevity  of  the  launch 
period,  negligible  liquid  spinup  is  expected  in  the  XM736  so 
that  the  angular  momentum  of  the  projectile  is  about  l\.%  less 
than  that  of  a comparable  solid  projectile  having  the  same 
total  mass  and  exterior  configuration.  Thus,  in  terms  of 
solid-projectile  behavior,  the  effective  axial  moment  of 
inertia  of  the  XM736  is  approximately  less  than  its  solid 
counterpart.  Pitch  inertia  is  negligibly  affected  by  liquid 
rotation  in  this  system,  as  previously  shown. 
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Figure  2.9.  Estimated  Spin  Versus  Time  for  the  Liquid  in  a Type  XM736  Projectile  at  Two  Values  of  Viscosity 


Exterior  Ballistic  Differential  Effects 


An  attempt  to  assess  the  exterior  ballistic  effects  of 
the  indicated  change  in  rotational  inertia  relative  to  a 
"ballistically  matched”  solid  projectile  was  made  in  the 
following  manner.  A set  of  runs  was  made  with  a modified 
point-mass  program  [/jJ  in  which  the  inertial  characteristics 
of  a projectile  having  characteristics  similar  to  the  M106 
projectile,  and  termed  the  standard  eight-inch  projectile, 
were  changed  incrementally  as  shown  below  in  Table  2.1. 

Shifts  in  range  and  deflection  relative  to  those  of  the 
standard  projectile  are  noted.  Projectile  characteristics 
are  shown  in  Table  2. 4* 

The  combined  effect  of  movement  of  the  center  of  gravity 
toward  the  nose  by  0.12  inch  and  reduction  in  the  axial 
moment  of  inertia  by  \$>  is  to  change  the  deflection  by 
about  1.36  milliradians.  This  magnitude  is  less  than  one 
deflection  probable  error  for  the  M106  system  [2,$]. 

Change  in  range  is  negligible.  Similar  ballistic  sensiti- 
vities for  the  M106  and  M509  projectiles  are  shown  in 
Tables  2.2,  2.3A,  and  2.3B. 

[2]  Op.  Cit. 

[I4J  Lieske,  R.P.,  and  Reiter,  M.L.  Equations  of  Motion  for 
a Modified  Point  Mass  Trajectory,  BRL  Report  No.  T^UTT 
Ballistics  Research  Labs,  Aberdeen,  Md. , March  1966. 

[53  Firing  Tables  for  Cannon,  8-Inch  Howitzer,  M2,  M2A1 
and  MI4.7  Firing  Projectile,  HE,  M106,  FT  8-J-3,  Hdqts 
Dept  of  Army,  October  I960. 


TABLE  2.1.  BALLISTIC  SENSITIVITY  OF  A STANDARD  PROJECTILE 
IN  M2  8-INCH  HOWITZER  AT  MAXIMUM  RANGE  TO  CHANGES 
IN  CG  POSITION  AND  AXIAL  MOMENT  OF  INERTIA 


Loo . CG 
(cal) 

Axial  MI 
( slug  ft* 

Rang© 

)(a) 

A Rang© 
(a) 

Def  1. 
(m) 

A D©f 1. 
(m) 

2.5 

0.5 

16,943 

0 

665 

0 

2.485 

o.5 

16,943 

0 

672 

7 

2.5 

0.14.8 

16,945 

2 

634 

- 31 

2.485 

0.14.8 

16,945 

2 

642 

- 23 

TABLE  2. 
IN  M2 

2.  BALLISTIC  SENSITIVITY  OF  THE  M106  PROJECTILE 
HOWITZER  TO  A JOINT  CHANGE  IN  CG  POSITION  AND 
AXIAL  MOMENT  OF  INERTIA 

QE  45° 

Vo 

Loc.  CG 

Axial  MI 

Rang© 

Defl. 

A Defl. 

(f/s) 

(cal) 

(kg  m3) 

(a) 

(a) 

(a) 

1380 

2.840 

0.553 

11751 

210 

0 

(Z5) 

2.825 

0.530 

11752 

201 

- 9 

1950 

2.840 

0.553 

16795 

417 

0 

(Z7) 

2.825 

0.530 

16796 

400 

- 17 
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TABLE  2.3A.  BALLISTIC  SENSITIVITY  OP  THE  M509  PROJECTILE 
IN  M2  AND  M110E2  HOWITZERS  TO  A JOINT  CHANGE  IN  CG 
POSITION  AND  AXIAL  MOMENT  OP  INERTIA 

Qg  14.5° 


V 

(f/s) 

Loc . CG 
(cal) 

Axial  MI 
(kg  m2) 

Range 

(m) 

Def  1 . 
(m) 

A Defl. 
(m) 

M2  Howitzer 

1906 

3.  £92 

0.570 

170£9 

298 

0 

( 25) 

3.  £77 

0.£47 

170£9 

285 

- 13 

M110E2 

Howitzer 

1040 

3.  £92 

0.570 

8 £47 

136 

0 

(Z3) 

3.  £77 

o.£47 

8 £48 

130 

- 6 

i960 

3.  £92 

0.570 

17622 

398 

0 

(Z7) 

3.  £77 

o.£47 

17622 

382 

- 16 

2440 

3.  £92 

0.570 

228£3 

621 

0 

(Z9) 

3.  £77 

o.£47 

228£3 

£9£ 

- 26 
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TABLE  2.3B.  BALLISTIC  SENSITIVITY  OP  THE  M509  PROJECTILE 
IN  THE  M110E2  HOWITZER  TO  A JOINT  CHANGE  IN  CG 
POSITION  AND  AXIAL  MOMENT  OF  INERTIA 


QE 

V /Zone 
o 

Modif.  Range 

Def  1 . 

A Defl. 

Correction 

(deg) 

(f/s) 

No 

."  (m) 

(m) 

(m) 

( mils ) 

12 

ioij.o/z3 

0 

3822 

11 

0 

0 

1 

3822 

11 

0 

0.1 

2k 

0 

6627 

k2 

0 

0 

1 

6627 

kl 

- 1 

0.3 

k$ 

0 

Q5k7 

136 

0 

0 

1 

8$kB 

130 

- 6 

0.7 

12 

1960/Z7 

0 

9399 

55 

0 

0 

1 

9399 

53 

- 2 

0.2 

2k 

0 

13967 

158 

0 

0 

1 

13966 

152 

- 6 

O.k 

kS 

0 

17622 

398 

0 

0 

1 

17622 

382 

- 16 

1.0 

12 

2^0/Z9 

0 

12869 

8ii 

0 

0 

1 

12669 

81 

- 3 

0.2 

2k 

0 

18369 

2k7 

0 

0 

1 

18368 

237 

- 10 

0.5 

k$ 

0 

22853 

621 

0 

0 

1 

22853 

595 

- 26 

1.2 

For 

modification 

number  n: 

n 

Loc.  of  CG 
(cal  re  nose) 

Axial  MI 
(kg  m2) 

0 

3-592 

0.570 

1 

3-577 

0.5k7 
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TABLE  2.4A.  CHARACTERISTICS  OP  THE  STANDARD  PROJECTILE 


caliber 

203.2  mm 

mass 

200  lb 

eg  position 

2.5  cal  aft  of  nose 

axial  moment  of  inertia 

0.5  slug  ft2  (0.6779  kg 

m2) 

pitch  moment  of  inertia 

4.0  slug  ft2  (5.4234  km 

m2  ) 

center  of  pressure 

1.73  cal  at  Mach  1.77 

projectile  length 

4.3  cal  (34*4  in) 

muzzle  velocity 

1950  f/s  at  zone  7 in  M2 

! howitzer 

initial  spin 

735  rad  sec-'*' 

TABLE  2.4B.  DRAG  COEFFICIENT  FOR  THE  STANDARD  PROJECTILE 


Mach  No. 

CD  (form) 

CD  (skin) 

Cp  (total) 

0.0 

0.126 

0.056 

0.182 

0.8 

0.126 

0.049 

0.175 

0.9 

0.190 

0.048 

0. 238 

1.0 

0.302 

0.047 

0.349 

1.1 

0.307 

0.046 

0.353 

1.2 

0.300 

0.046 

0.346 

1.5 

0.  262 

0.045 

0.307 

2.0 

0.210 

0.043 

0.253 
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TABLE  2.5A.  CHARACTERISTICS  OP  THE  M106  PROJECTILE 
FIRED  PROM  THE  M2A2  CANNON 


caliber 

203 • 2 mm 

mass 

200  lb 

eg  position  re  nose 

2. 8I4.O  cal 

projectile  length 

k-375  cal 

axial  moment  of  inertia 

0.553  kg  nT 

pitch  moment  of  inertia 

U-.2 70  kg  m" 

muzzle  velocity 

1950  f/s  at 

zone  7 

1380  f/s  at 

zone  5 

initial  spin 

735  rad/s  at 

zone  7 

520  rad/s  at 

zone  5 

TABLE  2.5B.  DRAG  COEFFICIENT*  FOR 

THE  M106  PROJECTILE 

Mach  No. 

CD 

0.00 

0.125 

0.75 

0.125 

0.85 

0.129 

0.90 

O.ltt-O 

0.95 

0.152 

1.00 

0.351 

1.05 

O.tt-OO 

1.10 

0.ij.00 

1.50 

0.356 

2.00 

0.305 

2.50 

0.280 

BRL  estimate  [ 6] 

[6]  Dubin,  J.A.,  et  al.  Ballistic  Similitude:  8 Inch  Ammu- 
nition, (SECRET),  Technical  Report  ij.165,  Picatinny 
Arsenal,  Dover,  N.J.,  June,  1973- 
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TABLE  2.6A.  CHARACTERISTICS  OP  THE  M509  ICM  PROJECTILE 
USED  WITH  THE  M2A2  AND  THE  XM201  CANNONS 


caliber 

203 • 2 mm 

mass 

20^.9  lb 

eg  position  re  nose 

3.592  cal 

projectile  length 

5- 674  cal 

axial  moment  of  inertia 

0.570  kg  m:; 

pitch  moment  of  inertia 

I4..  7676  kg  ms 

muzzle  velocity  (max) 

1906  f/s  in  M2A2 
22k 0 f/s  in  XM201 

initial  spin 

718.5  rad/s  in  M2A2 
1150  rad/s  in  XM201 

TABLE  2.6b.  DRAG  COEFFICIENT*  FOR  THE  *#09  PROJECTILE 


Mach  No. 

CD 

0.00 

0.130 

0.75 

0.130 

0.85 

O.II4.O 

0.90 

0.155 

1.00 

0.300 

1.05 

0.360 

1.10 

0.360 

1.50 

0.317 

2.00 

0.27  k 

2.50 

0.239 

BRL  estimate  [6] 

[6]  Dubin,  J.A.,  ©t  al.  Op.  Cit. 
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CHAPTER  III 

VIBRATION  OP  THE  LIQUID  IN  A 
SPINNING  LIQUID-PILLED  PROJECTILE 

Introduction 


In  an  effort  to  assess  the  consequences  of  vibration 
of  the  liquid  surface  while  in  Configuration  B on  the  flight 
stability  of  the  projectile,  a simple  methodology  is  pro- 
posed. First,  one  must  examine  the  admissible  shapes  which 
the  liquid  surface  can  assume.  For  a set  of  vibrational 
modes  of  interest  one  can  then  estimate  the  associated 
natural  vibrational  frequencies  for  the  liquid,  treated  as 
a conservative  system.  To  be  analytically  tractable  this 
treatment  will  assume  the  liquid  to  be  rotating  at  a constant 
angular  velocity  w.  Actually,  of  course,  the  liquid  does 
not  have  a constant  w everywhere  during  spinup  and,  further, 
at  any  point  in  the  liquid  w depends  upon  time.  Therefore, 
the  degree  of  credibility  of  results  derived  from  the 
above  assumption  will  depend  upon  the  relative  magnitude 
of  liquid  acceleration  due  to  spinup  and  the  centrifugal 
acceleration  due  to  spin.  If  the  latter  is  much  larger 
than  the  former,  it  is  plausible  to  treat  liquid  vibra- 
tion pseudostatically . 

The  ultimate  goal  of  the  present  analysis  is  to  com- 
pare the  natural  vibrational  frequencies  of  the  liquid  with 
the  frequencies  of  precession  and  nutation  of  the  entire 
projectile.  If  any  of  the  vibrational  frequencies  were 
found  to  remain  close  to  the  precessional  or  nutational 
frequency  of  the  projectile  during  flight,  a resonant  con- 
dition could  occur  in  which  the  system  vibrational  modes  excited 
each  other  at  their  common  frequency.  Conceivably  this 
could  cause  projectile  flight  instability  if  the  liquid 
vibration  was  severe  enough.  At  the  very  least  a "mode 
lock"  of  this  sort  would  increase  projectile  dispersion. 


-45 


To  motivate  further  developments,  one  should  note  that 
the  processional  and  nutational  frequencies  of  concern  are 
rather  low,  lying  in  the  band  from  0 to  20  hertz,  approxi- 
mately. For  stable  projectiles  such  as  the  M106  or  the 
m£09,  the  precessional  frequency  is  typically  about  1 to  3 
hertz  thruout  flight.  The  Appendix  includes  derivations  of 
the  equations  for  precessional  frequency  and  nutational  fre- 
quency of  a spin-stabilized  projectile. 

By  contrast  to  the  very  low  precessional  frequency, 
nutational  or  yawing  frequency  is  of  greater  concern  rela- 
tive to  projectile  stability“in  liquid-filled  projectiles 
since  this  frequency  is  such  that  a liquid  vibrational 
frequency  will  cross  it  during  spinup.  To  illustrate  the 
band  of  nutational  frequencies,  Figure  3.1  displays  the  nu- 
tational frequency  versus  time,  for  several  firing  zones, 
during  the  flight  of  the  M£09  projectile  from  the  XM201 
cannon.  Figure  3*2  shows  a similar  result  for  the  M106 
projectile  from  the  M2A2  cannon. 

A stable  projectile  is  one  in  which  pitching  or  yawing 
motions  are  ultimately  reduced  in  amplitude  during  flight 
without  the  projectile  being  completely  overturned.  In 
practice  a condition  of  neutral,  dynamic  stability,  such 
that  amplitudes  remain  constant,  is  difficult  to  obtain. 
If  the  system  starts  to  progressively  increase  its  yaw, 
it  does  so  quickly  and  ultimately  overturns. 
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Figure  3.1.  Nutational  Frequency  of  the  M{?09  Projectile  at  Several  Zones  in  the 
XM201  Cannon 


4b 


Figure  3.2.  Nutational  Frequency  and  Stability  Factor  During  Flight  of  the  M106 
Projectile  from  the  M2A2  Cannon 


Vibrational  Modes 


We  proceed  with  the  analysis  outlined  above  by  con- 
sidering perturbations  to  the  liquid  surface  of  Configuration 
B.  In  Figure  3.3  below,  the  equilibrium  position  of  the  free 
surface  of  the  liquid  at  its  inner  radius  is  designated  y . 
The  perturbation  or  deflection  of  this  surface  is  § where 


6 = f(x) 


(3.1) 


Thus  the  position  of  the  surface,  y,  is  given  by 


y = yG  + 6 . 


(3.2) 


1 


1 yo  | yU) 


i-  -s> 


^ * 


0 


*m 

T 


Figure  3*3.  Surface  of  the  Vibrating  Liquid 


The  liquid  volume  v is  given  by 


(3.3) 
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(3.14-) 


v = ff 

Then  with  a constant  liquid  volume  given  by  (3.I4.)  and  with 
5 small,  i.e., 


6 « y. 


/ 


V2 

-x  / 2 

171 


f(x)  dx  at  0 


(3.5) 


Equation  (3.5)  places  a constraint  on  the  form  of  f (x) . 
Now  expand  f(x)  in  a Fourier  series,  i.e., 

00  CD 

^ / \ \ ' 2 rr  n x \ ^ , 2 rr  n x t ^ \ 

f(x)  ■ ) a^  cos  — — + ) bn  sin  . (3.6) 

n«l  n=l 


Furthermore  since  we  are  interested  in  vibrational  modes 
which  might  amplify  yaw  and  destabilize  the  projectile,  con- 
sider only  odd  components  of  f(x)  for  the  present.  In  this 
special  case 

f(x)  . [ bn  sin  -L-n-  - . (3.7) 

n-1 


Applying  the  constraint  given  by  (3*5) > 


x /2 


— x / 2 

in 


^ , 2 rr  n x ^ 

b sin  ? r 0 

n l» 


n=l 


After  exchanging  the  order  of  operations 
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or 


~ r x / 2 
Y b f sin 

^ - n J -xj2 


2 if  n x 

— Z 


n=l 


0 


2E' 

nsl 


rr  n x 

b cos  — 7 — £ _ o 
n L 


(3.8) 


for  all  integer  n.  For  this  expression  to  hold  for  an  arbi 
trary  set  bn,  the  argument  of  the  cosine  must  be  an  odd 
integer  multiple  of  tt/2. 

it  n x 

- s (2  k • 1)  j , l^k 

, 1 < n 


Therefore, 

2 n x 

L = 2L 

2 k - 1 


and  with  (3.7) 


f(x) 


sin 


it  (2  k - 1)  x 
*m 


(3.9) 


(3.10) 


Since  we  are  principally  interested  in  low  frequency 
vibrations,  we  restrict  the  following  development  to  the 
fundamental  mode,  i.e.,  to  k = 1.  Then 

f(x)  = b sin  < x < ~ . (3.11) 
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The  Fourier  amplitude  b will  serve  as  the  single  degree  of 
freedom  characterizing  the  vibrations  of  the  liquid  surface. 


Energy  Considerations 

At  this  point  a brief  excursion  in  the  development  is 
necessary  to  discuss  the  pressure  in,  and  energy  of,  the 
rotating  liquid. 

The  rotational  kinetic  energy  of  the  liquid  is 
r x«/2  ~ 7m  t 

T , - tr  p / dx  f ra  r j u"  dr  (3.12) 

-V2  J y 

where  liquid  density  is  p and  where  the  angular  velocity  oo 
may  depend  upon  both  x and  r.  At  dynamic  equilibrium 

7 z y0.  , 

oo  r and 

o 

Trot  « | p “o  ym  (1  - q4)  <3'13) 

with 

q = y0/ym  • 


The  potential  energy  of  the  liquid  in  an  equilibrium 
configuration  is  due  only  to  the  compression  of  the  liquid 
under  centrifugal  force.  The  compression  of  the  differen- 
tial volume  element  dv, 

dv  - 2 it  r dr  dx  , (3*14) 
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is  due  to  the  pressure  needed  to  support  the  liquid  below 
r,  i.e.,  at  a radius  smaller  than  r.  In  a liquid  rotating 
as  a solid  body  with  constant  angular  velocity  03 


(3.15) 


which  upon  integration  yields 

<d*  P 

p(r)  s — ~ (r3  - y3)  , yQ  « r $ ym  . (3.16) 


If  k is  the  volumetric  compressibility  of  the  liquid,  i.e., 

k - i a v 

“ v a p * 

then  the  compressional  potential  energy  associated  with  the 
differential  volume  dv  is 

d Vc  - | p3  dv  . (3.17) 

Typically  k is  10“^  cm3  dyne-^  for  liquids. 

The  total  compressional  potential  energy 

Vc  = /T  I dv  • 

And,  (3.1U,  3.17)  yield 

V B it  k f d x C p3  r dr  (3.18) 

c J-x  /2  Jy 

m J 
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For  a liquid  in  dynamic  equilibrium,  p(r)  is  given  by 
(3*16).  Using  this  result 


To  appreciate  the  significance  of  the  value  of  liquid 
pressure,  the  rotational  kinetic  energy,  and  the  compressive 
potential  energy,  we  present  the  following  numerical  example. 

Example 

Using  values  used  in  previous  examples, 
w0  = 9^2  rad  sec”^ 

p - 1 gm  cm 


o 


or 


with 


’ - Vym 


(3.19) 


x s 55.88  cm  (22  in) 

TT1 


ym  s 7*62  cm  (3  in) 


yo  = 2.516  cm  (0.991  in) 


and 
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k ■ 10“ 10  cm3  dyne 
q = 0.3302  . 


Prom  (3*16),  the  pressure  at  radial  position  ym  is  given 


by 


2 pL  2 

co-  p y_ 

P(ym)  « 2 (1  - q8) 


m (1  > 0. 33023) 

2 


p(y  ) = 2.295  107  dyne  cm"3 
m 

- 22.65  atm 
= 332.95  psi 
And  from  (3.19)* 

vc  = 2T  (9^2)2  (10~10)  (55.88)  (7.62) 6 (1  - 0. 3302s)3 

V - 7.973  107  dyne  cm  or  7-973  joules  . 
c 

Finally,  from  (3*13) 

Trot  * £ <942)“  (5S>88)  <?-62>4  (1  - 0.3302^) 

T . a 1.297*4-  1011  dyne  cm 
rot 

T . * 12,9714-  joules  or  9,5&9  ft  lb- 
rot  x 
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The  potential  energy  of  the  liquid  at  equilibrium  is  only 
0.061$  of  the  rotational  kinetic  energy  under  this  condition. 
Therefore  one  would  not  expect  the  exchange  of  energy  between 
kinetic  and  compressional  potential  forms  to  contribute  sig- 
nificantly to  liquid  vibrations.  Compressional  potential 
energy  is  assumed  negligible  in  subsequent  calculations. 

Estimate  of  a Fundamental  Vibrational  Frequency 

At  this  point  we  return  to  the  principal  arguments 
associated  with  the  derivation  of  an  expression  for  the  vi- 
brational frequency  of  the  liquid  surface.  To  develop  this 
expression,  an  estimate  of  the  vibrational  kinetic  energy 
of  the  liquid  will  be  required. 

An  estimate  of  the  kinetic  energy  associated  with  a 
longitudinal  vibrational  mode  of  the  liquid  can  be  made 
simply  by  making  these  assumptions: 

(1)  The  liquid  surface  during  vibration  remains  axi- 
symmetric;  i.e.,  circumferential  modes  are  not  excited. 

(2)  The  functional  form  describing  the  liquid  surface 
changes  with  time  only  thru  time-dependence  of  the  coeffi- 
cients in  a Fourier  transform  of  the  function,  i.e.,  only 
thru  the  Fourier  amplitudes. 

(3)  The  column  length  for  flow  of  a liquid  element  at 
position  x (relative  to  the  center  of  the  disturbance)  is 
proportional  to  x. 

By  the  first  and  second  assumptions  the  first  odd  vi- 
brational mode  — 

y s b sin 

m 

has  time  derivative 
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(3.20) 


y - b sin 


IT  X 


At  the  surface  element 


2 rr  70  &*■  > 


the  path  or  column  length  of  a control  volume  within  which 
this  surface  element  vibrates  is 


c x 


0 < x < 


(3-21) 


by  the  third  assumption,  where  c is  a constant.  This  con- 
stant can  be  evaluated  by  requiring  that  the  volume  integral 
of  all  elemental  control  volumes  produces  the  volume  of  the 
liquid  v . That  is 

v 

vl  = / 2 nr  yc  c X ta 


or 


4 v, 


rr  y xr 
J o m 


(3.22) 


But 


v a rr  y*  x , 
l o 


(3.23) 


where  xQ  is  the  length  of  the  liquid  cylinder  in  Configuration  A, 
Then 


C 3 


1*0*1 
xm  y0 


(3.21;) 
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In  previous  examples 


x - 19.6  in 
0 

x s 22  in 

XU 

70  z 0.991  in 

ym  = 3 in 

so  that 

c - 1.471  • 


The  kinetic  energy  of  vibration  of  an  elemental  volume  is 
tf  c p yQ  x dx  y2  , 0 < x ^ 

The  total  vibrational  kinetic  energy  is  obtained  by  volume 
integration  of  this  expression  with  y given  by  (3.20). 

x /2 

Tvib  = * c P yo  b3  / m x sin*  (3.25) 

^ o m 

Tvib « 0 » y0  4 (3-26) 

16  rr 


or 


Tvib  = Kt  b* 


with 

Kt  » 


(TT3  + 4) 
16  TT 


c p y0 


(3.27) 
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Thus  the  vibrational  kinetic  energy  is  simply  proportional  to 
the  square  of  the  time  derivative  of  the  Fourier  amplitude. 

An  expression  for  the  vibrational  potential  energy 
associated  with  this  mode  is  developed  as  follows.  Let 
p(r,yo)  be  the  pressure  in  the  liquid  at  radial  position  r 
when  the  liquid  is  in  Configuration  B.  An  expression  for 
p(r,yQ)  is  given  in  (3.16).  Then,  providing  the  displacement 
6 of  th©  surface  from  yQ  is  small,  the  work  performed  to 
effect  a change  in  liquid  configuration  from  the  equilibrium 
position  yQ  to  a terminal  position  y is 

x / 2 b 

W = / m d x / 2 it  y [p  (y  + 5,  y - 5)  d6] 

J XmO  J Z — O 


with 


5 r z sin 


IT  X 


m 


(3.28) 


IT 


W = 5 p "S  *1  b 


(3.29) 


But  the  vibrational  potential  energy  is  equal  to  the 
work  done  to  change  the  configuration.  Thus 


vib 


= K b' 


with 


*o  *m 


(3.30) 


For  b equal  to  1 cm  and  the  other  parameters  having  the 
values  given  in  the  previous  example,  equals  I4.9 * 3 joules 

or  only  0.38$  of  the  rotational  kinetic  energy. 
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Having  obtained  expressions  for  the  potential  and 
kinetic  energy  of  vibration  of  the  liquid  in  terms  of  the 
Fourier  amplitude  b and  its  first  derivative,  one  can  obtain 
the  equation  of  motion  of  the  liquid  surface  by  direct  appli- 
cation of  Lagrange's  equation  for  a conservative  system. 


d 

as 


a t 


vib 


a b 


a v, 


vib 


a b 


r 0 


(3.3D 


This  result  does,  of  course,  neglect  the  effects  of  dissipa- 
tive forces. 

But  from  (3*27) 


a t 


vib 


a b 


= 2 Kt  b 


and  from  (3.30) 


a v. 


vib 


ft  - 2Kv 


Therefore , 

K 


b + gr—  b — 0 * 


(3.32) 


The  undamped  angular  frequency  associated  with  this 
vibrational  mode  by  inspection  of  (3*32)  is 


®vib  = < W 


(3.33) 
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^vib 

wo 


■[ 


(3.34) 


(ff*  + 4)  c xj 


Sine©  the  natural  vibratory  frequency  is  , 


The  value  of  «0  in  this  expression  is  interpreted  as  the 
effective  angular  velocity  of  the  liquid,  i.e.,  that  uniform 
angular  velocity  which  produces  the  observed  angular  momentum 
when  w is  not  uniform  thruout.  At  the  muzzle  spin  for  zone  7 
of  the  M2  howitzer 

uQ  a 735  rad  sec”1 


(3.35) 


With  the  parameter  values  used  in  previous  examples 


yo  s 0.991  in 


x - 22  in 

Cl 


C a 1.471 


« 0.06645  «0  * 


(3.36) 


and 


^y£b  3 hz  • 


For  zone  7 in  the  M110E2  howitzer 
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and 


w0  s 923  rad  sec 


-1 


Z 6l.3  hZ  • 

Using  (3.36)  and  the  time -dependent , effective  angular 
velocity  for  the  XM736  projectile  shown  in  Figure  2.9,  the 
value  of  the  function  V t)  for  this  system  has  been  com- 
puted. This  result  is  displayed  in  Figure  3*^*  Also  shown 
here  for  comparison  is  the  nutational  frequency  at  zone  7 
for  the  XM736  projectile. 

It  is  noted  that  the  spinup  of  a liquid  with  constant 
kinematic  viscosity  of  10  centistokes  occurs  rapidly  enough 
to  cause  the  vibrational  frequency  to  quickly  cross  over 
the  nutational  frequency.  In  this  case  liquid-vibration- 
induced  projectile  instability  is  unlikely.  With  a liquid 
having  lower  viscosity  while  in  Configuration  B,  the  cross - 
over  of  frequencies  is  not  so  abrupt  and  may  at  least  pro- 
duce additional  dispersion.  Considering  the  magnitude  of 
the  perturbing  effect  of  liquid  vibration,  complete  pro- 
jectile instability  is  unlikely. 
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Figure  3.4.  Comparison  of  Oscillatory  Modes  of  the  XM736  Liquid-Filled  Projectile 
(Zone  7 at  QE  45°  in  M110E2  Howitzer) 


In  view  of  the  approximations  required  to  perform 
this  type  of  analysis,  stronger  conclusions  are  not 
warranted."  An  experimental  program  to  better  define  the 
values  of  the  pertinent  parameters  characterizing  the  liquid 
certainly  is  suggested.  Additionally,  a program  to  define 
the  ballistic  dispersion  of  projectiles  with  low-viscosity 
liquid  fills  is  indicated. 


*-Note  [3.1]: 

A more  elaborate  analysis  such  as  a finite-element 
description  of  the  liquid,  embedded  in  a six-degree-of- 
freedom  flight  simulation  appears  to  be  a profitable 
direction  to  take  analytically;  however,  this  method  will 
likely  blur  some  of  the  insights  afforded  by  the  simpler 
procedure  used  here. 

Some  of  the  restrictive  assumptions  used  in  our 
analysis  could  be  removed.  As  an  example  we  note  that 
our  analysis  in  Chapter  III  neglects  certain  forces  arising 
because  of  the  non-Newtonian  character  of  the  coordinate 
system.  A coordinate  system  which  rotates  with  the  pro- 
jectile is  assumed  here.  Since  this  is  a non-Newtonian 
frame,  a modified  form  of  Euler's  equations  must  be  used. 
This  formulation  introduces  the  usual  centrifugal  and 
Coriolis  forces  as  well  as  an  "angular  acceleration" 
force  which  arises  from  the  angular  acceleration  of  the 
projectile.  The  latter  force  is  in  the  same  direction  as 
the  Coriolis  force  and  is  proportional  to  the  angular 
acceleration  of  the  projectile. 

For  a liquid  element  of  mass  m in  this  rotating  co- 
ordinate system,  the  apparent  acceleration  vector  a in 

— m 

the  moving  frame  due  to  an  external  force  F is  given  by 
the  equation 
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Not©  [3*1]  (continued) 


0)  x r_ 

— -Htna 


- (x) 


XV 

-m 


03  X 


(to  x P^) 


is  the  origin  of  the  moving  frame  relative  to  a 
Newtonian  frame 

is  the  position  of  the  mass  element  in  the  moving 
frame 

is  the  apparent  velocity  of  the  mass  element  in 
the  moving  frame 

is  the  angular  velocity  of  the  moving  frame  in  a 
Newtonian  frame 

is  the  angular  acceleration  of  the  moving  frame 
(projectile) 

The  term  - go  x (w  x r ) is  the  centrifugal  acceleration 
- - -m 

which  has  been  treated.  The  term  - 2 w x p is  the  Coriolis 

— — ma 

acceleration.  In  a right-handed  frame  with  the  x-axis  along 
the  spin  axis  of  the  projectile,  this  term  reduces  to 

• 2 a (j  k • i j) 

with  and  k unit  vectors  in  the  y-  and  z-directions . If 
the  transverse  components  of  velocity  are  small,  this  term 
is  negligible  and,  in  fact,  was  neglected.  The  term 
- a)  x r is  the  angular  acceleration  component  of  accelera- 

— —Cl 

tion  and  reduces  to 


where 


So 


r 

-*n 


Sma 


to 


0) 


-kwy+iwz 

This  term  has  also  been  neglected  due  to  the  small  value  of  <0. 
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APPENDIX 

DERIVATIONS  AND  COMPUTER  PROGRAMS 

Estimate  of  Processional  and  Nutational  Frequencies 
of  a Spin-Stabilized  Projectile 

For  a gyroscopic  ally  stable  spinning  system  with  angu- 
lar velocity  o)  and  axial  moment  of  inertia  I , the  pro- 
cessional angular  velocity  9 obtained  when  the  gyro  feels 
a torque  t is  given  by 


9 


(Al) 


If  the  torque  is  produced  by  aerodynamic  forces  where  the 

center  of  pressure  is  - XBrn  calibers  in  front  of  the  eg,  then 

sm  0 


t r - X D C„  A q s - X D CM  •<  A q 
sm  N sm  N 


( A2a) 


where,  notationally  and  by  definition, 

D m caliber 

Cjj  s normal  force  coefficient 
C„  a normal  force  derivative  coefficient 
•<  * projectile  angle  of  attack 
A m reference  area 


q ■ dynamic  pressure 


(A2b) 


(A2c) 
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Sine©  the  overturning-moment  coefficient 


CM  = - X_  CM 
M sm  N 


<Xog-Xcp)  CN 


equation  (A2a)  may  be  written 


(A3) 


T S CM  D A q 

T = CM  ■<  D A q (A4) 

The  natural  precessional  frequency  at  constant  -<  is 


V S 


S? 

CM  D A 

3TT 


co 


(A5) 


Yaw  Frequency 

The  solution  of  the  linearized  equations  for  the  yawing 
motion  of  a projectile  is  given  by  McShane  et  al  on  p.  6i|4 
of  [7)»  The  solution  involves  superposition  of  two  sinusoids 
having  different  frequencies.  The  two  frequencies  are 

21,2  = (Ito)  (rad/s) 


or 

^1,2  * 2tT  21,2 


(hz) 


( A6) 


with 

o = (l-s"1)^ 


(A7) 


where  s is  the  gyroaeoplo  stability  factor  given  by 

s - JA 

‘ 8 IB  y D3  , 


[7]  MeShane,  E.,  Kelly,  J.  and  Reno,  P.  Exterior  Ballistics, 
University  of  Denver  Press,  c.  1953* 
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or 


s z 


0) 


Tl 


XB  CM 


D3  q 


(AS) 


The  frequencies  V ^ ^ have  been  called  the  "nutational" 
and  "processional"  frequency,  respectively.  However,  it 
should  be  noted  that  is  not  the  reciprocal  of  the  time 
between  yaw  maxima.  Because  of  the  superposition  of  the 
high-  and  low-frequency  components  of  motion,  the  frequency 
with  which  maximum  yaw  occurs  is  the  difference  frequency: 


(A9) 


This  is  called  the  yaw  or  nutational  frequency  in  this  report. 
An  alternative  derivation  for  V starts  with  an  expression 
for  the  wavelength  of  yaw  A,  in  calibers.  Then, 


Vy  = BTf  (h.) 


(A10) 


with  projectile  velocity  V and  caliber  D. 

An  expression  for  A is  obtained  from  p.  65l  of  [7], 
Using  our  notation 


A = 


2 ; tb  v 

IA  D co  o 


Then,  (A10)  and  (All)  yield 


*1  A 

Vv:m 


B 


(All) 


(A12) 


This  frequency  has  been  computed  as  a function  of  time  for 
several  combinations  of  projectiles  and  cannons.  Results 
for  the  M509  projectile  in  the  XM201  cannon  are  shown  in 
Figure  3 • 1 • Comparable  results  for  the  M106  projectile  in 
the  M2A2  cannon  are  shown  in  Figure  3-2. 

[?]  McShane,  E.,  Kelly,  J.  and  Reno,  F.  Op.  Cit. 
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Derivation  of  an  Equation 
for  Spin  Decay  in  Projectiles 


Glossary  of  Terms 

N z projectile  spin  (rad/sec) 

N a initial  spin  (rad/sec) 
o 

V z projectile  velocity  (m/sec) 

V = muzzle  velocity  (m/sec) 
o 

t - time  since  launch  (sec) 
p s air  density  (kg/m  ) 

D a projectile  caliber  (m) 

M s projectile  mass  (kg) 

p 

IA  z projectile  axial  moment  of  inertia  (kg  m ) 
= spin  damping  moment  coefficient 
Kp  z zero-lift  drag  coefficient 
k s velocity  decay  parameter  (m”'*') 
x - rangewise  coordinate  (m) 


A 

For  flat  trajectories  at  low  QE,  we  assume  negligible 
gravitational  effects. 


(A13) 


2 


V 


2 


V 


M 


(AII4.) 


(A15) 


with 

Then, 


with 


d(v)-1  KD  P D 

“TT" = M 


v"1  - v"1 

O 


KD  p D‘ 

— H — 


V . (V"1  + kt)"1 


k - 


kd  P d 

— H — 


(A16) 

(A17) 

(Alb) 


Ns  - V N 


with 


P D*  K. 


-c  s 


N 

IT 


V"1  + kt 
o 


d(ln  N)  = - dt 


V-*  + kt 
o 


In  N 


N 


N 


. - £ in  (v;1  + kt) 


t 

0 


i N 
ln*o" 


i in  (1  + k Vo  t) 


N , 

r»u^v  o t)! 


^ ka  M D 

p = i - 

K *D  A 


(A19) 

(A20) 

(A21) 

(A22) 


(A23) 

(A24) 

(A25) 
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Vi  it  - s- 

(i  + k vo  tr 

Integration  of  (A17)  produces 

V = V e-ta 
o 

And  with  (Al6), 

t = k"1  (V-1  - V”1) 

or 

t r k-1  (V”1  ek*  - V”1) 
o o 

t = (e1^  - 1)  / (k  Vo) 

Substitution  of  (A29)  into  (A26)  yields 

N = N e-'3  kx  . 
o 


(A26) 


(A27) 


(A28) 


(A29) 


(A30) 


An  Example 


For  long  trajectories  the  variation  of  drag  coefficient 
and  spin-damping  moment  coefficient  with  Mach  number  renders 
the  above  results  quite  approximate.  However,  in  some  in- 
stances this  approximation  may  be  adequate.  In  this  example 
analytic  results  for  spin  damping  versus  time  are  compared 
with  those  produced  by  a computer  simulation  in  which  varia- 
tion with  Mach  number  is  considered. 

Take  the  M106,  8 inch  HE  projectile  as  fired  from  the 
M2  howitzer.  The  maximum  range  trajectory  will  be  considered. 
In  this  case  an  average  altitude  ASL  is  about  10,000  ft. 

At  this  altitude  air  density  is  about  0.7^  sea  level  stan- 
dard. Thus 

p = 0. 74.  PQ  = 0.9065  kg  m-3 
Po  = 1.225  kg  m-3 
Other  parameter  values  are 

N(0)  = Nq  = 735  rad  sec-1 
Vq  a 59ll-.il-  m sec-1 
C„  = 0.30  (effective)  or 
Kd  = 0.1178  (effective) 

= 0.006  rad-1  (effective) 

D = 0.203  m 
M = 90.72  kg 
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Then 


I = O.678  kg  m" 

k = kd  P D2  M'1 
k = ij..8506  10~^  m-1 

and 

“ = P »k  ka  II1 

_£  -1 

x = 1.3622  10  ^m 
(3  = -</k 
3 = 0.2808 

V k = 0.02883  sec-1 
o 

Prom  (A26) 

N = Nq  (1  + k Vq  t)_t3 

N = 735  (1  + 0.02883  t)"0,28°8  (A31) 

The  result  in  (A31)  is  plotted  in  the  following  graph 

with  selected  variables  from  a simulated  trajectory.  For 
some  purposes  the  agreement  shown  between  the  analytical  es- 
timate and  the  more  exact  simulation  may  be  satisfactory. 
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Estimate  of  Stability  of  Spin-Stabilized  Projectiles 
Having  Oscillating  Inertial  Properties 


The  objective  of  this  section  is  to  display  and  illus- 
trate by  examples  a numerical  technique  for  evaluating  the 
stability  of  spin-stabilized  projectiles  which  are  subject 
to  oscillations  in  the  position  of  the  center  of  gravity 
(X  ) and  simultaneous  oscillations  in  the  transverse  moment 
of  inertia  (I  ).  To  preserve  the  greatest  generality  of 

J 

approach,  the  equations  of  motion  describing  projectile 
pitch-yaw  dynamics  are  solved  numerically"  in  the  time  do- 
main with  eg  position  and  transverse  inertia  given  by  the 
following  functions: 


Xcg  * Xcgo  + 6Xcg  sln[2’,t  (*o+*'t)  + 

Iy  s Iyo  Sly  si.nl  2xt  » ( A32) 


with  X 


ego* 


6Xcg’  Jyo'  6V 


V constants  for  t £ 0. 


The  solution  is  carried  far  enough  in  time  to  determine 
whether  the  nutational  amplitude  is  decreasing.  While  ad- 
mittedly a brute-force  approach  such  as  this  lacks  the 
elegance  of  a frequency-domain  approach  to  stability,  there 
is  no  need  to  linearize  or  to  assume  the  amplitudes  fiXcg 
end  5ly  are  small  or  that  the  inertial  properties  oscillate 
at  a constant  driving  frequency.  Further,  some  insight  is 
gained  regarding  the  behavior  of  the  projectile  as  the 
driving  frequency  matches  the  nutational  frequency  of  the 
projectile . 


* Using  a fourth-order  Runge-Kutta  procedure  with  time 
step  0.002  sec. 
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Using  principally  the  notation  of  AMCP  706-165 , reference 
[8],  the  following  auxilia^  variables  in  the  equations  of 
motion  are  defined. 


H - ^T5T  [CN^  - 2 CD  - (CMq+V] 


Pa  A d ,_-2 


M 3 - V - - k CM  m.  a A q d CM  I 
2 m y y 


-1 


q * 2 pa  v 


p Ad 
ft  . rc 
TmT  CN 


cd  + k;2  cmp  1 


p . ix  i;1  w a 


a 


d 


P g 


a v 


-1 


m d 


m d 


(A33) 


In  this  notation: 

A ■ reference  area  of  the  projectile 
d m caliber  of  the  projectile 
m ■ mass  of  the  projectile 

[8]  Engineering  Design  Handbook:  Liquid-Pilled  Projectile 
Design,  AMCfr  706-165,  April 
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I = longitudinal  moment  of  inertia 
I • transverse  moment  of  inertia 
k s longitudinal  radius  of  gyration  in  calibers 

A* 

k - transverse  radius  of  gyration  in  calibers 
7 

v a velocity  of  the  projectile 
a)  m spin  of  the  projectile 
density 

component  of  gravity  normal  to  v 
zero-lift  drag  coefficient 
normal  force  derivative  coefficient 
overturning  moment  derivative  coefficient 
a the  pitch  damping  moment  coefficient 
Magnus  moment  derivative  coefficient 

Using  a right-handed  coordinate  system  as  in  [8]  with 
the  x-axis  coinciding  with  the  form  axis  of  the  projectile, 
positive  toward  the  nose,  the  pitching  angular  motion  in  a 
vertical  plane  (about  the  horizontal  transverse  axis)  will 
be  denoted  by  and  the  yaw  about  the  mutually  orthogonal 
axis  will  be  denoted  by  \jr.  With  this  notation,  the  equa- 
tions of  motion  are: 

[8]  Op.  Cit. 


pn  s air 

oL 

g s the 
-L 

CD  ■ the 

CN  s the 

CM  s the 

Cvr  + 

Mq  \ 

CM  b the 
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J a a-2  M t)  - a”2  P T \V  - a”1  H J - a”1  P ^ + a”2  G 
\|/  • a”2  P T t>  + a"2  M y + a”1  P J - a”1  H \|/  . (A34) 

Adopting  the  systematic  notation 


= * , (A35) 

equations  (A34)  become 

x s A x + b (A36) 

with 

X a [ X-j^  Xj  Xj^  ] 

b s [0  0 b3  0]  ' . 

The  elements  of  the  A matrix,  {a.,},  are  given  below 
-2  1 J 

in  units  of  sec 

all  a a12  = allj.  * 0 
a13  = 1 

a21  = a22  s a23  * 0 
&2k  ~ 1 
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a31  = 

a“2  M 8 A q d CM  i”1 

y 

a32  = 

- a"2  P T 

■ 

a33  = 

- w V-1  d2  4 A q i;1  [0  - CD  + k^2  C ) 

V 

- a"1  H 

S 

- A q v"1  m'1  [0  - 2 0 - k"2  (C  +C  )) 

V * q * 

a3/*  = 

- a”1  P s - I I"1  w 

^ j 

V a 

- a32 

a/j.2  = a31 


%3  = 

“ a3l|. 

f" 

II 

a33 

b3  * 

a“2  G . I I"1  g v”1  (A37) 

**  y jl 

Examples 

Using  the  equations  (A36)  and  (A37)«  two  projectiles 
have  been  treated  as  numerical  examples.  Reference  [8J  pro- 
vides data  on  the  solid,  WP  loaded,  1$2  ram  XMlj.10  projectile. 
Under  conditions  in  which  a portion  of  the  white  phosphorous 
fill  has  liquified,  both  the  effective  transverse  moment  of 
inertia  and  center  of  gravity  can  be  expected  to  oscillate 
in  flight.  Whereas  the  amplitude  of  these  oscillations  may 
be  slight,  a persistent  resonance  of  the  oscillations  with 
the  nutational  frequency  of  the  projectile  can  cause  the 
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nutational  amplitude  to  continuously  increase.  At  a launch 
Mach  number  of  1.5»  the  nominal  frequency  of  yaw  maxima  for 
this  system  is  17*5  hz.  Accordingly  several  numerical  ex- 
periments were  performed  in  which  VQ  (in  (A32))  was  set  to 
17.5  and  V was  set  to  0.0.  During  these  experiments  6X 

eg 

and  ol  were  varied  systematically  to  determine  the  region 

V 

of  stability.  The  absolute  stability  criterion  of  diminish- 
ing yaw  amplitude  was  used  here.  Preliminary  experiments 
indicated  that  stability  is  adversely  affected  when  §X 

eg 

end  5 1 are  in  phase.  Thus,  all  experiments  were  mm 

j 

under  the  worst-case  phase,  namely  for  * in  (A32)  set  to 
zero. 

A second  example  was  selected  for  comparitive  purposes. 
The  numerical  values  of  this  example  are  best  estimates  of 
the  parameters  of  the  XM73&  liquid-filled  projectile  at  a 
launch  Mach  number  of  unity.  Axial  moment  of  inertia  re- 
flects only  that  of  the  metal  parts,  indicating  that  the 
angular  momentum  of  the  liquid  at  launch  is  treated  as 
negligible.  The  amplitudes  of  the  inertial  increments  were 
selected,  somewhat  arbitrarily,  by  taking  values  proportional 
to  the  differences  observed  in  the  properties  of  Configura- 
tions A and  B of  Chapter  1.  The  values  of  the  parameters 
used  in  both  examples  are  displayed  in  Table  Al.  Proce- 
dures for  examining  the  stability  region  for  the  second 
example  were  identical  to  those  employed  for  the  first. 
Stability  was  examined  at  a forcing  frequency  equal  to  the 
nominal  nutational  frequency  of  7*77  hz.  Additionally, 
runs  were  made  in  which  the  frequency  was  swept  linearly 
from  6.0  hz,  at  the  rate  of  0.25  hz/sec,  to  8.5  hz  at  10  sec. 
These  runs  clearly  showed  that,  for  certain  values  of  qX 

eg 

and  the  projectile  will  remain  quite  stable  under  a 

condition  of  moving  forcing  frequency,  whereas  the  projectile 
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will  become  unstable  when  forced  at  a constant,  nutatlonal 
frequency.  The  stability  region  for  both  examples  is  dis- 
played in  Figure  A2.  Under  the  assumed  conditions,  the 
projectiles  for  both  examples  are  stable. 


j 
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TABLE  Al.  PARAMETER  VALUES  FOR 
GYROSCOPIC  STABILITY  ANALYSIS 


parameter 

symbol 

value 

XM410  XM736 

dimension 

caliber 

d 

0.5 

0.6667 

ft 

proj.  mass 

m 

1.313 

6.31 

slug 

long,  inertia 

Xx 

0.0446 

0.4036 

slug  ft2 

trans.  inertia 

*y 

0 . 154-8 

3 • 5164 

slug  ft2 

amplitude  of 
incr.  in  trans. 
inertia 

sly 

0.0006 

O.OI83 

slug  ft2 

nuzzle  velocity 

V 

o 

1675.5 

1117 

ft  s”1 

proj.  spin 

G) 

526.4 

573.4 

rad  s"1 

air  density 

pa 

2.3769 

10”3 

slug  ft”3 

drag  coef. 

CD 

0.50 

0.30 

normal  force 

CN^ 

2.90 

2.10 

rad"X 

pitch  damping 

-5.00 

-4.60 

(rad  sec”X)”X 

magnus  moment 

CM 

V 

0.30 

-0.10 

rad 

center  of  press. 

XcP 

1.40 

1.108 

cal 

proj.  center  of 
gravity 

Xcg 

1.85 

3.472 

cal 

amplitude  of 
incr.  in  c.g. 

fiXcg 

0.04 

0.020 

cal 

gravitational 

component 

gx 

0.0 

0.0 

ft  s”2 
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Figure  A2.  Stability  Regions  for  Projectiles  in  Examples  1 and  2. 


Computer  Program  for  Numerical  Solution 
of  One-Dimensional  Equations  for  Spinup  of  Liquid 

The  source  programs  shown  on  the  following  pages  were 
written  in  the  FORTRAN  IV  language  for  the  IBM  360-65  com- 
puter. Comments  in  the  listings  introduce  each  main  pro- 
gram and  subprogram  describing  its  function  and  delineating 
the  principal  operations  and  variables. 

Following  each  program  is  a sample  of  the  output  pro- 
duced by  the  program. 


Next  page  is  blank. 
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2 
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I 3 
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21 
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28 
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31 

32 
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34 
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36 

37 

38 

39 


$ . lOR  • MGEoRGE  • fKP  = 29tult>ILS  = S6U*  T lhE  = 3bO , PAGES  = S0  0 
C****  DYNAMICS  or  SPIN  up  of  LIQUID-FILLED  projectile 
c 

DIMENSION  TITLE*  (20)  fU(22)  **P<88)  tX(II) 

COMMON  AMR0A*N*X»DXtCl »C2t  TO 
READ  TITLE  ANq  j.MPUT  PARAMS 

READ  (StlOO)  TlTLEtDT»TO*TCOH?DX*XOf WL IMt NtNPRINT 
lOO  FORMAT (20a4/6FAO,0»2I2) 

WRITE  ( 6 1 SO  0 ) NtTCGN 

5o  0 FORMAT ( 1mo,3HN  *» I3»2X»6HTCON  = » F I 0 • 2) 

C**»*  PRINT  COLUMN  HEADINGS 

WRITE  (6*200)  TI  TLE*DT»DXtXO*WLIM. 

200  FORMAT (IHI,?0A4/1H0IIHTIME  STEP  =tEI5.5tI4H  SPACE  STEP  =*EI5,5* 

1 2X4HXO  =tFl5,5*2X6H^LlM  = * E I 5 • 5/ 1 HO * 9X * I H I * 6X4HU < 2 ) * 6X4HU ( 3 ) * 

2 ftX4nU(4)  » ft  X4HU ( 5 ) *6X4HU(6)  *6X4HU(7)  *6X4HU<8)  *6X4HU(9)  * 

3 SXSMU(IO) .SXSHU(II) ,5X*SH0MEGA* 10H  LlU  SPIN) 

C 

c«***  COMPUTE  constants 

NMisn-I 
X { I ) s XO 
SUM*  xO/2 • 

A M R D A 3 I • / 13  X * * 2 
00  2 J = 2*N 

X(J)sDX*FLOAT(J-l)+XO 
SUM=SUM+X( J) 

2 CONTINUE 

SUMs^UM  X ( »| ) / 2 • 

B=I ,/SUM 
C2*0, 0288*TCOM 
L'l  =-C2*0  • 2808 
C 

C***»  INITIALIZE  STATE  VECTOR 
U(  1)=WLIM 
00  4 vJ  = 2 » N m 1 
U(  J)  =k»LIM 
4 CONTINUE 
U(N) si .0 
C 

C**»*  INITIALIZE  RUNGE-KUTTA  SUBROUTINE 
C#***  SOLVE  DIFP  eons  FOR  DERIVATIVE  ICS 
T = 0 . 0 

CALL  DlFEO(TtUf4) 

19  K0UNTs V 
C 

C#***  START  OF  SOLUTION  LOOP 

20  CONTINUE 
C 

C**«*  MOVE  STATE  FROM  T TO  T*DT 

CALL  KUTTA(T.OT»U,^Pf Hf2f0lFE0) 

C 

C#***  COMPUTE  EFFECTIVE  ANGULAR  VELOCITY 
SUM=o.O 
00  30  J s I t Nf  1 1 

SUM=SUM* (U(J)*U(J*1) )*(X(J)*X(J*1) )/4# 

30  CONTINUE 

omega=b*sum 

FSPIn=92J.*0MEGA 
TP=TCON*T+TO 
IF  (TP,GT • ftft , ) CALL  EXIT 
KOUNT*KOUnT ♦ 1 
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40 

IF  (KOU  JT.cD.N»RI«*4T)  GO  TO 

4 1 

if(T.gt-  i.o)  call  exit 

42 

GO  TO  20 

43 

40 

WRITF  (O* JOG) TP* (U(I) *]=2*11) *0*tGA*r SPIN 

44 

300 

FORMAT ( 1 H * 13F1 0*5) 

45 

IF(T.GT.I.O)  CALL  EXIT 

46 

Go  TO  19 

47 

End 

48 

r 

SUh ROUTINE  OI FEU ( T I MF * U * KU T T A ) 

L 

differential  equations  for  angular  vel.  in  cone 

49 

Dimension  u(22)  *x(ii> 

50 

COMMON  AMHDA*N*X*DX*CI *C2* TO 

fi 

U ( N ♦ 1 ) = A MOD  A*  (U  (2)  -U  ( 1 ) )*(  (X  (2)  /X  ( 1 ) ) **2d  . ) /2. 

*2 

NM  .1  =N-  1 

53 

WO=l  . 

54 

uo  lo  I =2  * NM 1 

*5 

**1=  ( 1 ,-OX/X  ( I ) ) 

c6 

*P1=1 .♦OX/X ( I ) 

c7 

U ( I 4N)  sAMHDA*  (U  ( I-l ) *WMl-2.*tJ  ( j ) *W04U<  !♦!  j *WP1 ) 

58 

10 

CONTINUE 

*9 

U(  *N)=0.0 

6 0 

RETURN 

6 1 

ENU 

SUBROUTINE  KUTT A ( T * Dt * V t W t NEQ t NORDP 1 , D I FFEQ ) 

DIMENSION  V(NEQ*NDRDP1>  t W ( NEQ t NOR DPI « 4) 

DT2*DT*0#5 
DT6*oT/6.0 
DO  1 I*1»NEQ 
DO  1 J*1*NdRDP1 

1 W<I»J*1)«V(I*J) 

DO  2 K*  1 * J 
L*KM 

GD  TO  (3*3*4) « K 

3 DTW3DT2 
GD  TO  5 

4 DTW«DT 

5 TW«T*DTW 

DO  6 I»1*NEQ 
DO  6 J«2*N0RDPl 
Jl * J-l 

WPsW(I«Jltl) ♦ « (I* J, K)*DTW 
V(Ifjl)«WP 

6 W ( I * Jl *U>  *WP 

CALL  DIFFEQ(TW*V*K) 

DD  2 1*1*  NEQ 

2 W(I*mORDP1*L)*V(I*NOROP1) 

DO  7 J*2*NORDPl 

Jl*J-l 

DO  7 I*1*NEQ 

7 V ( I * Jl ) *W  ( I * Jl  * 1 ) ♦DT  6*  CWCIt  Jtl)*2'0*CH(l9  JtZ)  ♦VCX«J»3n*tf  (X  t J*4)  ) 

T*T  W 

CALL  D IFFEQ ( T * V * 4 ) 

return 

end 
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DYNAMICS  OF  SPINUP  of  liquid  CONFIGURATION  b*  case 
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Computer  Program  for  Modified 
Point-Mass  Exterior  Ballistics 


The  program  used  to  determine  the  ballistic  differential 
effects  of  changes  in  inertial  parameters  is  a general-purpose, 
exterior  ballistics  program  capable  of  simulating  the  flight 
of  gun-boosted  rockets  as  well  as  purely  ballistic  systems. 

The  program  employs  two  alternative  ways  of  treating  the 
aerodynamic  forces  normal  to  the  body.  The  first  option  simply 
assumes  trailing  behavior,  that  is,  that  yaw  is  always  zero. 
This  option  is  inadequate  for  the  purpose  of  this  study.  The 
second  option  computes  an  equilibrium  yaw,  or  "yaw  of  repose," 
necessary  to  process  the  velocity  vector  in  the  vertical  plane 
at  the  proper  rate.  Using  the  yaw  of  repose,  normal  body 
forces  are  computed  and  resolved  into  components  in  an 
inertial  frame  of  reference.  These  are  added  to  the  drag, 
gravitational,  and  Coriolis  forces  to  complete  a point-mass 
description.  This  procedure  gives  satisfactory  agreement 
with  experiment,  providing  accurate  aerodynamic  data  exist 
and  that  the  projectile  displays  adequate  stability.  This 
option  is  exercised  by  setting  switch  IOPTY  equal  to  unity. 
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»S I *02750440 


C 

~C~ 

C 

C~ 

c 

~c~ 


“EXTERIOR  BALLISTICS  OF  BOOSTED  ROCKETS 
A THREE-OEGREE-OF-FREEOOM  MOOEL  APPLICABLE  WHERE 
TRAILING  OR -FOLLOWING  BEHAVIDRCAN  BE  ASSUMED 


REAL  RS(40i) fTS(40n  *VS (4Ql ) »CJRT ( 1 1 ) iCJVT(ll) 

DIMENSION 'TrTLE(2Ol»UU2rrwPT4Br*XMT0L(II)  ,CDTBL  ( llTt COT0LU1T 
ItTMACH(ll)  #TKA(U)  .TKOYAW(ll)  • TKL  (11)  * TKM  ( 1 1 ) , TKF  (11)  • 

2 TKT(ll)  fTKH(lir»TKSTm  tTCPTTl)* 

INTEGER  *2  CHAR(l)/***/ 


00000100 
GOOD 02 00 
00000300 
00000400“ 
00000500 

Hjoowsiro — 

00000700 

oooooDOir  ~ 

00000900 

OOOOIDOO 

00001100 


ASSIGN  CONSTANTS  NEEOEO  BY  DIFFERENTIAL  EQUATIONS  TO  COMMON 


COMMON  EMOtEMBtSPI tFC»BR ATE fOELT#TO»TBi ISWtVi THETA tFFCTRtC ALSO f 
T“  »-SLTrR  tTENU>  CNACH  f REYNLITt  RES  I S » C AL  f DL  ONG  • lOPTYf  TAW  # AMON  f 

2 BMOMtPSlt WTAREAt ISEP t NABLEt IGLIOE 
CUM  NGN“ 


00001200  - 
00001300 

0 0*0  0T4W” 
00001500 
OO0OT60G 
00001700 
“OGDOiinnr- 


/SRCOM/TM 

COMMON/COFCOM/XCG»SMARGtEM, THIOtPRNU  * ALTRIM,CNATRM,GLI0E,EPSTHE  00001900 
T »5TAFAC#YWNU#TK^tTK0YAW»TKLlTKKTTKF»TKT»TKH>TKStTCPtTMACH#-NARTBLOOOO2(rOC 

common/wincom/rwc»xwc  00002100 

* ~xmtbltcdtbl*, cotbl  »ntbl 


COMmON/DRGCOM/ 

common/snocom/caoens 


00007200 

00002300 


OOOO2400-- 
00002500 
00002600 
00002700 
00 0028 (TO  * 
00002900 


C****  TABLES  OF  AEROOYNAMIC  COEFFICIENTS  IS  PARAM,  SET  INPUT  SET  -1. 
c»##*  rr  PARAMETERS  NTBT  AND  NARTBL  ARE  POTH  ZERO#  ENOOGENOUS 
C**#*  FUNCTIONAL  FITS  TO  THE  AEROOYNAMIC  T ABLES (WITH  THE  T387  FORM) 
C*»»*  WILL  BE  USED*  SEE  SUBROUTINE  ACOEFS.  

C****  IF  ONLY  NARTBL  IS  ZERO,  THE  ZERO-LIFT  ORAG  TABLE  IS  REQUIREO 


C*"***- WTTHOUT~REUUIR  ING  TABLES ‘ FOR  THE~ OTHER  AERO * COEFF I CTENTSi 00003000 

C#***  IF  CERTAIN  AERO  COEFFICIENTS  ARE  DEFINEO  (KNOWN),  THESE  00003100 

C»*##  CAN  PE  READ  WITH  THE  OTHERS  LEFT  BLANK*-'  THE*  PROGRAM  MILL 000832110 

C**#*  USE  THE  TABULATEO  COEFFICIENTS  ANO  DEFAULT  TO  THE  ENDOGENOUS  00003300 

FUNCTIONS  FOR  THOSE  ENTERED  AS  ZERO*”  ~ 00003400 


0*PR0JECT ILE  CALIBER,  MILLIMETERS 
"EMU*  INITIAL  PRDJECTTLtT'HAPS * LBH 


PARAMETER  INPUT  1 00003500 
INPUT  2^000 03  6TTCT 


~c~* 

c 

~TT 

C 

c 

~TT 

C 

IT 

C 

”C~ 

C 

~C 


EMB-BURNT  MASS,  LB* 
FCbNOHINAL “THRUST  LEVEL  ,*“LFF~ 


INPUT  3 00003700 
INPUT  4 O0003B00 

SPI*SPECIFIC  IMPULSE  OF  ROCKET  PROPELLANT,  LBF/LBM/SEC  INPUT  5 00003900 

— BRATEWPROPELLANT  BURNING  RATE!  LBM7SEC ENDOGENOUS  VARIABLE0 (TOO 4000 

OELT-THRUST  RISE  TIMEt  SEC  INPUT  6 00004100 

~ TDWIGNTTTON  TIME  FGR*  NDCKEY"MOTOR ENDOGENOUS  ' VAKI ABLE00O042OO 


IN  SUBROUTINE  »BURNf  THE  THRUST  OECAY  TIME  IS  ASSUMEO 

“A*  TYPICAL  VALUE  0 ^T~5EC* 


00004300 

00004400 


EQLUL  TO  THE  THRUST  PISE  TIMET 
TB*EFFECT I VE  BURNING  INTERVAL,  SEC  ENOOGENOUS  VAR  I ABLEO 0004500 

A SWITCH  SIGNALING  COMMENCEMENT  OF  BURNING  ENDG* 


ISWi 

IENO-A  SWITCH  SIGNALING  ENO  OF  BURNING 
V*PRGJECTTLE "VELOCITY , M/Sec 

VO-MUZZLE  VELOCITY  OF  THE  PROJECTILE,  FT/SEC. 
THETA* ATTITUDE  OF*PRO  JECTILE,  ~DEG  — 


ENOO*  VAR I ABLEO 00 04700 
END 0 . ’ V A H I A B LE  000 0’48'OTT 
INPUT  7 00004900 


“97- 


VARIABLE00005100 


CALSQ«CALIBER  SQUAREOf  H««2 ENOO, 

HJM.TELTTCITTTJF^REaUwIND",  "“W/SEC  REa13“TN~TW  TT/SFCT  ENDU* 

VWF«VEL0CITY  OF  HEAOWINU  In_FT/SEC  

H0« ‘INITIAL'  ALTITUDE,  FT* 
hterm-terminal  altituoe,  ft, 

FFcTR-FOrM  FACTOR~ReLATIVE  TO  PROGRAMMED  DRAG  FUNCTION, 

SEE  FUNCTION  *ORAG«  FOR  THE  SPECIFIC  ORaG  FUNCTION  USEO. 

GEO-TNTTTAL  UUa  DR  ANT"ELE  NATION**  DEg.  ” ~ ~ 

OQE«OUADRANT  ELEVATION  INCREMENT,  OEG*  INPUT  1300005900 

step-time  step  in  Numerical  integration  procEoure,  sEcV  input  T400006000 

TM-TIME  AT  WHICH  BURNING  OF  ROCKET  MOTOR  SHOULO  COMMENCE.  INPUT  1500006100 
NQE-NUMBER  OF  INCREMENTS  of  QUADRANT  ELEVATION  INPUT  1600006200 

NPRINT«NUMBER  OF  TIME  STEPS  EXECUTEO  BETWEEN  PRINTS  INPUT  1700006300 

~CCnTTnDE  


INPUT  8 00005300 
INPUT  8 00005W 
INPUT  1000005500 
INPUT  Il0o805600 
00005700 
"INPUT  T80088S808' 


ooooTm 

ENOO*  VARIABLE00006500 
EnDD*.  VARIABLE0000660C 
CONSTAKT00006700_ 
CONST*  ANTo  0 00660  0 
INPUT  1800006800 


c 

c 

c 

c 

c 

—£ 

c 

c 

c 

c ~ 
c 

c 

c 

c 

C ' 
c 

~ IT- 
c 
c 
c 
c 
c 

c 

c 

c 

c 

c 

c 


ALT* TRUE  ALTITUOE  ABOVE  SEa  LEVEL*  M 
R-RAOIUS  FROM  CENTER  OF  EARTH  TO  PROJECTILE  , M 

re«nominal  raoius  of  the  earth  at  the  equatop,  h 
Omega-angllaR  velocity  of  The  Earth*  RaD/5ec 

IOPTY-  A SWITCH  INDICATING  CHOICE  OF  Yaw  OPTION, 

ropTY.  I ppooucrs  computation  uf~yaw  of  repose  Ftnr&imtm  ppojsctoootqm 

IOPTY-  0 SIGNIFIES  A TRAILING  PROJECTILE  WITHOUT  SPIN,  FOR  00007100 

This  Option  The  following  inputs  are  unnecessary.  d80d?2o8 

SPINO-INITIAL  SPIN*  RAO/SEc  INPUT  1900007300 

XCG-POSITION  OF  center  OF  GRAVITY  aFT  OF  NOSE*  CALIBERS  INPUT  2000007480*~ 
XCP-POS I T I ON  OF  CENTER  OF  PRESSURE  AFT  OF  NOSE*  CAL  ENOO.  VAR  I ABLE0 0 00750 0 

-clong«projecttle  ltnsht^tn  calibers,^ Twur^rodwsoo 

AMOM-LONGITUOINAL  MOMENT  OF  INERTIA  OF  THE  PROJECTILE,  KG*M**2  00007700 

Input  22oooo78dd 

BMOM-TRANSVERSE  MOMENT  OF  INERTIA  of  The  PROJECTILE*  KG*M**2  00007900 

INPUT  2300008000 

VCW-VELOCITY  OF  cROSSWINO  FROM  RIGHT  LOOKING  OOWNRANGE (REAO  IN  FT/0000B100 

jfcPUy  74$  0 008280 

WTAREA-WETTEO  AREA  RATIO  USEO  IN  COMPUTING  SKIN  0000B300 

FRICTION  ORAG  INPUT  2500008400 

ISEP*A  SWITCH  INOICaTING  CHOICE  OF  SEPARATE  00008500 

COMPUTATION  OF  SkIN  FRICTION  DRAG.  INPUT  26000G8600 

- 1 IF  FRICTION  ORAG  IS  COMPUTEO  SEPARATELY  ANO  ADOEO  TO  FORM  ORAG0000B700 
= 5"  IF  "FRICTION  ORAG  IS  INCLUDED  I N'DPaTT  FUNCTION, 

SMARG-PROJECTILE  static  margin,  cal. 

PSI-ANGULAR  ORIENTATION  of  Yaw  VECTOR 
srng-slant  range  to  projectile  POSITION,  m 
VHXF-VELOCITY  OF  LAUNCHER  IN  RANGEWISE  DIRECTION, 
vhyf-velocity  of  launcher  in  VERTICAL  oirection. 

" CONTINUE 


RWC  IS  HEAOWINO  COEF, 

XWC  IS  CROSSWINO  COEF* 

PSI  May  BE  COMPUTEO  BY  REMOVING  «C* 

IN  SUBROUTINE  FLIGHT. 

CaOENS  IS  THE  CORRECTION  FaCTOR  FOR  AIR  DENSITY  RELATIVE  TO  STaNOA00009900 

oo  0 1 oooo 


S FROM  COMMENT  CAROS 


“88088*801) 
ENOO.  VARIABLE00008900 
ENOO,  VaRIaQLE00009000 
ENOO,  VARIABLE00009100 
(FT/SEC)  00009200 

(FT/SEC)  00009300 

O0  00840  O' 
00009500 
00009600 
00009700 
00009800 


external  FLIGHT  00010100 

EQUIVALENCE  (U(l) tX), (U(2) ,Y) , (U(3) tZl , (U(4) ,SPIN),  00010200 
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1 <U<5) ?X0) ? <U<6) ?Y0) ? <U<7>  ?Z0>  ? (U<8) ?SRINO) ? 

zr~nj  ( 9 r VXDOTVUJ  ( 1 0TVYDTT)  r<UTllT  ?Z5DT?TUTT2T  ?SDD'J ~ 

read  in  run  description?  constants  in  flight  equations  ano 
INITIAL  C0N0IT10NS* 

READ  <5?256?END*30)  TITLE?NT8l?NARTBL 
256  FORMAT (20A4/2I2) 

nrmm.EQ7«T'GO^ro  255 : 

REAO  (5  ?250  ) (XMTtJL(I)  ?1*1?NT8L) 

READ  ( 5?  250 ) (CDTBL  (l)~t  1*1  ?NTBL) 

READ  (5?250 ) (CDTBL ( 1 ) ? I»1 ?NTBL) 

250  FORMAT(8F10*0) 


00010300 

00010400 

00010500 

00010600“ 

00010700 
0 001  WTO 
00010900 

00OTT000 

0001 1100 

00011700 

00011300 

00011400 


WRITE  (6?252)  TITLE  00011500 

~ ”252  F0RWAT(1H1  ?20A471H0  ? I0H  mACRTNO  ?10H  “COEF  DRAG?T0N  DRAG  I NCR)  00011600 

00  253  1*1 tNTBL  00011700 

WRITE  (6? 254)  XMTBL (IT?  CDTbL i I > ? COTBLTIT  — 0O“01T800 

253  CONTINUE  00011900 

IF(NARTBL.EQ.O)  GO  00  755  “00012055 

REAO  ( 5 ? 250 ) <TMACH(I) ?I*lfNARTBL)  00012100 

; RTAff~~(59?5'0r  rTRATI  )TT*1TN ARTffO  00012200 


( TKDYAW ( 1 ) » I* 1 1 NARTBL ) 
(TKLfI) ?1«1?NaRTBL) 

< TKM ( I ) ?I*1?NaRTBL>  • 
READ  (5?  250)  ( T CPTiytT *Tt NOT TB C) “ 
REAO  (5?  250 ) (TKF(l) *I*1?NaRTBL) 


REAO  (5?250) 
READ  ( 5? 250) 
READ  (5?250) 


00012300 


000 12400 
00012500 
00017600“ 
00012700 


READ“<5»250)  ( TKT  n)TT*T?  N AR  T BET"  

READ  ( 5 ? 250 ) (TKH ( I ) 1 1*1 *NaRTBL> 

READ  (5? 250 ) < TKSTITt  TiT ? NARTBL) 

254  FORMAT ( 1H  ?3F10.4) 

WRITE  T6?Z72) 

272  FORMAT ( 1H0  ? 1 OH  MACH  N0?6X?2HKA?5X?5HKDYAW? 

T 8X  ?ZHKLT8X77HKK71  OF!  ~ TPT  CALI  

00  257  I*ltNARTBL 

WRITE~  (6? 270)  TKACH ( I)  tTKATTl  ?TKDYAW UT?TKL  (T)  ? TKM  ( I ) ?TCP  (IT 
270  FORMAT ( 1H  ?6F10,5> 

257  CONTINUE”  ~ 

255  CONTINUE 


"PKDV13T0NA1  CONSTANT  GLIUE-aKGLE^UaN  T5 
C***#  SWITCH  IGLIOE  MUST  BE  SET  TO  1 FOR  CONST • GLIDE  ANGLE  TRAJECTORY, 
READ  (5? 264)  IGLID5? ALTRlKfCNATRM  *5ETDE?EPSTHE 
264  FORMAT  ( 1 1 ?9X?  4F1 0 *0 ) 

TF  n GLI DE • NE • D GO  TtTI 

WRITE  ( 6 ? 268 ) 

268  TURHATnHOTlSH TRIM*  ANBUEA~N* 6X? 9HCNA  ( TRIKT? 

1 15H  GLIDE  ANGLE?  R?15H  TOLERENCE?  «) 


WRITE  (6?266)  ALTRIM  ? CNA  TRM  ? GL IDE  * EP3THE 
266  FORMAT  (1H  ?4F15,5) 

PRDVTS1DNAL  CONSTANT  GCTD5  ANGLE  JAN  T5“ 

1 READ  (5?  2?ENO®30)  Tl TLE ?0 ?EMO ?EMB ? FC ? SP1 , OELT ? VO ? VWF ?HO? HTERM ? 

~ ~ 0~TFCTR?QEO?DQE?STEP?TW?NQE'.NPRINT?TOPTY 

2 FORMAT (20A4/BF10,0/7F10*0?3I3) 

SWITCH  NABLE  IS  SET  FROM  0 TO  ~TT  AT  TIME  TENABL 


00012800 

00012900 

00015000 

00013100 

00013200“ 

00013300 

00013400 

00013500 

00015600 

00013700 

00013800“ 

00013900 

000T4000 

00014100 

00014200 

00014300 

00014400 

00014500 

00014600 

00014700 
”00014500 
00014900 
”0  0 015000 
00015100 

00015200 

00015300 
00015400  ~ 
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258  READ  (5*260)  C AOENS , VCW * TENA0L » THI D 00015500 

2 hJT  FORMAT  ( 4F 1 0 . (T) DTOT^TTO 

WHITE  (6*262)  TEN  AB(_  * THI  D * CADENS  00015700 

262  FORMAT(lHO*l4HENABLE  TIME  • *T10.T*5X*  OO016dO6 

1 21HTHRUST  DRAG  FACTOR  » • FI  0 .4  * 5X * 1 8HA IR  DENS  FACTOR  ■ *F10.4>  00015900 

' IF(10PTY.N£.1J  GO  TO  25  OOOlStfOO 

RE AO  26*SPIN0*XCG*CL0NQ* AMOM*BMOM*WTAHEA* ISEP  00016100 

rETaO  "26  * VflXr*VHYF **wC*  XtfC 60016260 

26  FORMAT (6F10.0* 12)  00016300 

GO  fD  27  00016400 

25  SP IND*0 • 0 00016500 

XCG-0.  5o0l66*G 

CLDNG**0.  00016700 

~aMOM»0.  * ’ 00016800 

8MDM«0.  00016900 

' wTarEa*0 • OOTITOOO 

ISEP-0  00017100 

VHXF»(T.d  0OO1T21T0 

VHYF»0 • 0 00017300 

RWOO’.O ” ~ ’ 60D1T4DT5 

XWC«0.0  00017500 

SMARGbO.O  0001*600 

STaFAC«0.0  00017700 

Ya^NU-0.0  ToGlYBOO 

PRNU-0.0  00017900 

— PSr-OYfl  00013000 

OMY-O.O  00018100 

"C  START  QUaD-ELEV  LOOP  O0OIb2oO 

27  QE»QEO-OQE  0001B300 

SDO«0.0  0001B400 

CAL*0*1 .E-3  0001 8500 

DD  3_TQF«1*NQE  — 00 OTB 600 

QE»QE^DQE  00010700 

THE T A»QE /57 • 29578  06018300 

TD«1.E10  00018900 

EM«EMO  00  0 1 9007) 

C IEND  IS  A SWITCH  SIGNALLING  END  OF  BURNING  00019100 

TENO-0 " (TOO  192 00 

C ASSIGN  TIME  INCREMENT  FOR  INTEGRATION.  00019300 

c CO 01*400 

DT»STEP  00019500 

C NaBLE  IS  A SWITCH  SIGNALING  CONTROLS  DEPLDYEO  FDR  GUIDED  FLIGHT.  OO0196O0" 

NABLE-0  00019700 

T PRINT  anD^aSeL" "RUN  DT5GtfTFTIWr~C0NSTAHTS*  aW  ~ ; 0WI*9T0 

C INITIAL  CONDITIONS  00019900 

print  9*title*ffctr*vd*emo,emb*o*oe*ot,fc*spi »vwf*hd*hterm*vcw  T0626G66 

1 *RWC*XWC  00020100 

9 FDRMaT(1h120A4/1H010X5HFFCTR13X2HV013X2HM013X2HMB14X1hD/  OO0202Od 

1 1H  5F15.6/1H0*  00020300 

2"6X9HQUAD  ELEV8X7HTH  STEP9)T6HTHRUST5X1TJHSP  IMP(JL5F9X6HV^¥TND/  00020400 

3 1H  5F15.6/1H04X1 1HINIT  ALT  * FT4X 1 1 HTERM  ALT  * FT  1 5H  VEL  XWIND*FT/S*  00020500 

4 12X*3HRWC*12X*3HXWC*/1 H 5F15.6)  - 00020600 
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IEBUPDTr  LOG  PAGE  6065 


new  Faster 


PRINT  91 1 VHXF  t VHYF 

~FTJRHATT/0N  VHXF“i~TFr272t25H  FT75EC 
l/) 


VRYF 


00020700 


00020900 


92 

PRINT  92 tOELT t TH 

FORMAT ( 1 H 19HTHRUST  RISE  TIME  » tFl2.4#5H  SECt8H  TM  * 

00021000 
tF 12 • 4 t 00021  100 

C 

1 5H  SEC) 

00021200 

00021300 

z 

c 

yo»initial  altitude,  h 

DO  02 1 40  0 
00021500 

c 

c 

HO-ALTITUOE  READ  IN  IN  FT 

yterm»terminal  altitude*  m 

00021500 

00021700 

c 

HTERM- TERMINAL  ALTITUDE  READ  IN  IN  FT 

D0021500 

IMIOPTY.NE.1)  GO  To  29 

00021900 

28 

PfflNT  ZB,XC6,in.URUVAMUM,¥M0tf 

FORMAT (lHOfllHLOC  OF  CG  ■ *E1 0. * * 2X ,3HCaL, UH  PROJ  LENGTH 

00022000 
■*E10*4*  00022100 

"1  2X, 3HCAL,  15H  181  At  M~OF~T  i.EI1.5,2X*THKG  M**2, 
2 15H  TRANS  M OF  I -.Ell .5*2X,7HKG  M**2) 

00022360 

00022300 

29 

1TT5PI jEQ.0.0)  GO  TO  60 
BRATE-FC/SPI 

50022A58 

00022500 

TT( rRATETE Q.O. 01  GO  TO  60 
TB» ( EMO-EM0) /BRATE 

“66022600 

00022700 

40 

PRINT  AO » TB 

FORMAT ( 1H0 t 22HEFFECT I VE  BURN  TIME  » tF10.4*4H  SEC) 

(HID  2Z50  0 
00022900 

GO~TO  6F~ 

— 58023000 

60 

TB*0  • 

00023100 

C 

0RATF-O. 

* 60  0232  00 

00023300 

C 

61 

COMPUTE  AUXILLIARY  CONSTANTS  AND  REDIWENSI ON  INPUTS 
CALSQ*0**2*I.E-6 

06023500 

00023500 

ULONG»D»l . OE-5»CLONG 
VELO»0.3048*VO 

00023600 

00023700 

C 

SUPVEL  TS  THE  SUPREMIUM  "OF“PRUJECTILE  VELOCITY. 
SUPVEL»VELO 

00023900 

SUP ALT  IS" 5UPREHIUW  OF  PROJECTILE  ALTITUOE. 
SUPALTaO.O 

~ 00024000 
00024100 

VW»D«3048#VWF 

00024200 

VCW«0,3046*VCW 

00024300 

YO«OV3G4B*HO  ~ ' 
YTERM«0.3048*HTERM 

00024400 

00024500 

RTERHaRE* YTERM 

00024500 

c 

INITIALIZE  TIME*  X*  X-DOT*  Y AND  Y-OOT. 

00024700 

— c 

T»0. 

00024000 

00024900 

~X»0.  ~ ‘ 

0 0 0250  0 0 

Y«YO 

00025100 

z»o. 

SPIN-SPINO 

60025200 

00025300 

VHX« *5048* VHXF 

06025405 

VHY».3048*VHYF 

00025500 

YO»VELO*SIN ( THETA) ♦VHY 
THETA *57. 29578* A TAN  prD/XDr 


00025700 
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V-SORT (XD##2^YD##2 ) 

ZD-5. 

SPIND-0.0 
XDD-0 • 0 
YDD-0.0 
ZDD-0 . 0 
YAW-0.0 

aLt»yO  ~ 

SRNQ-0.0 

ISW-IBURN  (ft  ALT# THETA  f V) 
IF(ISW.EQ.l)  GO  TO  7O 
GO  TO  71 


00025900 

00026000 

00026100 

00026200 

00026300 

00026400” 

00026500 

50026600 

00026700 

0T626F0T~ 

00026900 

56~0275W~ 


70  TO-O.O  00027100 

DT-$TEP~74.  * 00027200 

00027300 

INITIALIZE  RUNGE-KUTTA  SUBROUTINE'  ' ~ 00027400 

71  CALL  RUNGE1 (U#WPt4f2t FLIGHT)  00027500 

SOLVE  FLIGHT ‘EQUATION?  FOR" INITIAL  CONDITIONS.  00027600 

00027700 

' CALL” FLIGHT <ttUt4)“'  00027800 

00027900 

INITIALIZE  COUNTER  FOR  DETERMINING  NUMBER  OF  POINTS  fO  BE  PLOTTED  00028000 

AT  END  OF  TRAJECTORY  SOLUTION.  00028100 

00fr28?00 

NPLOT-0  00028300 

" 7- ootedW 

00028500 

INITIALIZE  COUNTER  FOR  COUnTINGT^nES  PER  PaGE.  00028600 

, 00028700 

LINE«0  ^ “ 665288OT 

00028900 

1 P RTNT -•NPffI  NT  015W900TT 

YAWDEG-YAW«57.3  00029100 

GO  TO  PRINT  OUT  INITIAL  CONDITIONS^  00029250 

00029300 

GO  TO  4 00029406 


00029500 


START  or  SOLUTION- XOOF7  5aVE  CAST  VALUES  OE  'FLIGHT 

\ 

"VAN  I ABLEST 

‘ 06029600 

00029700 

5 XP-X 
RP-R 
ZP-Z 
XDP-XD 

, 

06029800 

00029900 

OOTOO0O0 

00030100 

* YDP-TD  ' 

» 

00030200 

VP-V 

00030300 

thetap-theta 

OOOSO^OO 

CMACHP»CMACH 

00030500 

RESlSP-RESIS 

000  3O6JT0 

YAWP»YAWDEQ 

00030700 

SPlNP-SPfN  ' ' ' 

' 00030800 

SMAR6P-SMAR0 

00030900 

STAFAP-STAFAC  0003IOTTO 
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NEW  MASTER 


I EBUPDTE  LOG  PAGE  0007 


PSIP-PSI 

SRNGP-SRNG  

call  runge-kutta  SUBROUTINE  to  solve  flight  equations  from 
T TO  T*OT. 

CALL  RUNGE2 ( T tOT ) 


C 

C 

C 

~C 

C 

"C 

c 

c 


SAVE  positional  COOROINaTLS  of  projectile  for  later  plotting  of 
trajectory. 

save  range  and  flight  time  in  arrays  for  subsequent 

POLYNOMIAL  fit 

SRNG*SQRT(X*X*  <Y-YO)**2*Z*Z) 

IFU.lt. TENABL)  go  to  520 
IF(NABLE.EQ.I)  GO  To  520 
“*  PRINT  530  * T 

530  FORMAT(1HO*18HENABLEMENT  TIME  * tFl0.3) 

NABLE*1 

IF  (NPLOT.EQ.400)  GO  TO  52o 
IF  (SRNG.GT.SRNGEM)  GO  TO  520 

nplot-nplotu 

TS(NPLOT)*T 
RS(NPLoT>*SRNG 
VS ( NPLOT) *v 

520  CONTINUE 

YAW0EG»YAW*57.3 
rF(V.GT.SUPVEL)  SUPVEL-V 
IF(ALT.GT.SUPALT)  SUPALT-ALT 
BB“IF(ISW.EQ.O)  GO  TO  50 
IF(IENO.EQ.I)  GO  TO  51 

IFU.GE.TO*TB*SELTJ  SO'TCTTS™ ~ 

GO  TO  51 

49  DT-STEP 
IENO«l 

CALL  BURN (T  * XMASS  * THRUST) 

PRINT  80.  XMASS.THRUST.ALT.VtT 

~B0  FORMATUH(T.9H  HASS  * •F10.47TIH'  ^THRUST  «~iF10V4. 

1 13H  ALTITUOE  « .F10.2.10H  SPEEO  « *F10.2.l0H  TIME  « .F10.3) 
GO  TO  51 

50  ISW«IBURN(T.ALT.THETA.V) 

IF (EMO.EQtEMB)  ISW«0 
IF(ISW.EQ.O)  GO  TO  51 


00031100 
DO  03 1 20D 
00031300 
00031400 
00031500 
00031000 
00031700 

oo'OMfioO- 

00031900 
00032000 
00032100 
00032200 
00032300 
00  0 324710 
00032500 
00  03260  0 
00032700 
00032000 
00032900 
00033000 
00033100 
00033200 
**00033300 
00033400 
00033500 
00033600 
00033700 
00033800 
00033900 
**00034000 
00034100 
0 0 0 3420  0 
00034300 
00034400 
00034500 
00034600 
00034700 
0(ro34fforu  — 
00034900 
“00035000 
00035100 
DO  03S200 
0 0 03530  0 
” W03547HJ 


E0 


TO«T 

PRINT  20.  TO 
FORMAT  UH0.15HBURN 


STARTS  AT  .F10.4.4H  SEC) 


00035500 
00  03SS00 
00035700 
~ 00TJ 35  87)0“ 
00035900 
~7J<nJ36U7Hr“ 
00036100 
“ 0UU362O7J 
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5T 

53 

4 

• *** 


DT-STEP/4. 

TP«r«fsrppnifTvi 

IFJIPRINT.EQ.O)  GO  TO  53 
“GO  tD  52 
IPRINT--NPRINT 


AOVANCE  LINES  COUNTER  ANO  CHECK  IF  TIME  TO  EJECT  PAGE  AND  LABEL. 


00036300 
60036460 
00036500 
001*366  ~00 
00036700 
~ 00036^50" 
00036900 

60037000“ 


00037100 

00037200 

00037300 


0 0037  ♦ oir 

00037500 


LINEfLINEM 

dmy«yawnu 

OHY  IS  a DUMMY  VARIABLE  USEO  FDR  OUTPUT  DF  CHOICE 
IF  (LINE.LE.O)  GO  To  6 
LlNE*-50 

“PRINTS  .TITLE  " " 00037600 

FDRMAT(lHl20A4/lH0f9HTIME»SECSt7X3HX#M5X5HALTtMt7Xt3HZtM2JU 00037700 

1 8HXDDT.M/$2X8HYDDT.M/S5X5HVtM/$lX8HMACH  NO . 2X7H0R AG • LB . 00037800 

2 9H  THETAtDt9H  YAW.0.9H  SPIN.R/S.9H  STa  FAC  »8H  OUMMY  V)  00037900 

PRINT  8 1 T .X.ALY.ZtXO.YD.V.  CMACH . RES I S »THETa . YAWDEG • SPIN. ST  aPAC  .OMYO 00380 00 
FORMAT (1H  lF9.3.3FlOrlt3F10.itF9.2t2F9#l,4F9.2)  00038100 

!^(T7GTV3bbV)  gdTTO  3<T 


52 


RULE  FOR  STOPPING  jpj-UTION  - STOP  WHEN  PRDJECTILE  HITS  GROUNO. 

IF(. NOT. (R.LE.RTERM. AND. THETA. LT. 0.0) ) GO  TO  5 
INTERPOLATE  SOLUTION  VARIABLES  FOR  R-RTERM 


00038*00 

00038300 


YE*YTERM  ~~ 

TE»T-DT*(R-RTERM)/(R-RP) 
DEL«<T-TE)/DT 
XE*X-DEL«(X-XP) 

ZE«2-DEL* ( Z-ZP ) 

XOE*XO»OEL# ( XO-XOP ) 

YDF-YD-OEL*  < YD-YDR!” 
VE«V-OEL<MV-VP) 

THeTaE«ThETa-OEL* ( THET A-THeTaP) 

CMACHE«CMACH-OEL*(CMACH-CMaCHP) 

RESISE-RESIS-DELMReSIS-RESISP) 

YAWE«YaWOEG-DEL«(YAwDEG-YAWP) 

SFTNE«SPI  N-DE1*T5RTN-SPTNP1 

STAFAE«STAFAC-OEL#  tSTAFAC-STAFAF 

SMARGE«SMARG-DEL*ISMARG-SMaRGP) 

PSIE«PSI-DEL«(PSI-PSIP) 

SRNGE«SRNG-DEL»(SRNG-SRNGP) 

XP LOT(NPLOT)bXE 
"TRLOTTl  •NPTOTT «TE~ 


00638400 

00038500 

~0W38660~ 

00038700 


00038800 

00038900 

W039966 

00039100 

06039200" 

00039300 


06039400 

00039500 


00 039600 

00039700 

06039800 

00039900 

00040000 

) 

00040100 

PRINT  DUT  SOLUTION  VARIABLES  FOR  Y.TTERM, 

PRINT  8.  TE.  XE.YE.TE.XDE.YDE .VE  t CMA£HE»TfESTSE*THETAE.Y4VE.~ 
1 SPINE. STAFAEiSRNGE 


SUP VEL-SUPVEL/O. 3048 
RANGE-SQRT (*E»»2*ZE*«Z1 


00040300 
90941140  0“ 
00040500 

00040600' 

00040700 
0004980  0~ 
00040900 
~ 9994198IT 
00041100 

— 09047290” 

00041300 
TTO047  409 
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RANGE “RANGE* ( 1 • ♦ ( RANGE/RE ) **2/6* ) / 1 852 . 

PRTFTTHA'XIWM  VETO C 1 TY *~R7VN GE iTANIT^U TlTUDEi 

PRINT  90  * SUPVEL*  RANGEtSUPALT 
90  FORMAT  (1H0 1 20HMAX~PROJ  VELOCITY  » *F15. 4 1 4HT7St 

1 3X«12HMAX  RANGE  * *F15.4*11H  NAUT  MI LES * 3X * 1 OHMAX  ALT  ■ 

2 F15.4,0H  METERS) 


~PLOT ‘TME'TRA JECT0RY~JU5T  COMPUTE 
LABEL  PLOT  WITH  TITLE*  QE  AND  VO. 


PRINT  10#TITLE*QEtVO 

~T0~FORMAT ~( IHO 15X20 A4/4H  QE*F5TIVI0H"_DEG*~ 
3 CONTINUE 


VCJVF6.  I * 4H~F75T~ 


~Z~ 

C 


CALL  POLFIT  <RS*TS*NPLOTt 3*0 *CJRT *6*6* .TRUE.* SOEV* 
120HFLIGHT  TIME  VS  RANGE) 

CALL  POLFIT (RS * VS  * NPLOT ♦ 3 * 0 *C JVT ♦ H*HINV * . TRUE . *SOEVt 

'20MPR0 J7  VEU7  TST  RANGE) 


00041500 

TTOGATSCTO 

00041700 
00041000 
00041900 
"00042000 
00042100 
W04220D 
00042300 
00042400  " 
00042500 
00042600 
00042700 

"W0T2W0 

00042900 

*00043000 

00043100 

00043200 

00043300 

OffiT4340D 


RETURN  FOR  ANOTHER  CASE.  00043500 

••*•«••«•*««#• ###•••»#»**«•**«»»«*»**•«•••••«••#••••••«#*##**#•»»• 000 436 oO 

GO  TO  1 00043700 

30  “CALL  CTI T — 00043000 

ENO  00043900 

■ SUBROUTINE  SOUND  ( A*G*RHUWTSGO>'  ' ~ ' 00044000 

COMMON  EMO*EMb*SPl tFCt0RATEtDELTtTOtTB»ISWtVtTHETA*FFCTR*CALSQ*  00044100 

1 VW  * VCW>  ALT  * R*IEND  *CMACW*  REYNLD t RESTS • CAL*  OLONG* I OPTY  * YAW*  AMOM*  00044200 

2 BMOMtPSI * WTAREA* I SEP  * NA0LE t IGLIDE 
C0MM0N/SN0C0M/CADEN5 

EQUIVALENCE  ( Y t ALT) 


TT 

C 

XT 

C 

C“ 

C 


SUBROUTINE  COMPUTES  THE  SPEEO  OF  SOUND  IN  M/SEC 
VERSUS  ALTITUDE  IN  HE  TER  ST  ALSO  COMPUTED  IS  THE 
ACCELERATION  OUE  TO  GRAVITY  IN  M/SEC/SEC  AND  THE 
AIR  DENSITY  IN  KG/M**T  AND'  THE  ABSOLUTE  VISCOSITY 
OF  THE  AIR  IN  KG/M/SEC.  NOTE  THAT  REYNOLD’S  NUMBER 


~GS9V026* ( 6 . 370E6/C6,TT0E6*YT 
0*6.356766E6*Y 
IFTY.LT.  II 01Bt  07}  -BO~TO  1" 
IF(Y.LE.20063.12)  GO  TO  2 
TFCY.LE.3ZI6I.9)~G0  TO~3  A 
IF(Y.LE#47350.(N>  GO  TO  4 
“TFrr.LEt5242BYBBT~GO  TO  5 
IF<Y.LE. 61591. 03)  GO  TO  6 
iFTTt LTrT9994VT4J  BO  TT TT 
RHO»0.4636*EXP(-0.12207E-3*Y) 

t»tempepature  in  degrees  kelvtn — 


00044300 
— 00044400 
00044500 

00044*00“ 

00044700 
~ DOO44B0O~ 
00044900 
0004500D~ 
00045100 
00045200 
00045300 
~ 7)0045400“ 
00045500 
00045600 
00045700 

00045H00 

00045900 
00046000  “ 
00046100 

00046200” 

00046300 

00046400 

00046500 

00045500’ 


TM0O.65 

~B~A*Ztr.  05T»5QRT  (TT 
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RHO«RHD*CAOENS 


— VT5C0-OYOO4“67*  (T* 1 1 0 . > * ( T/217.78) *#1 .b 
THIS  IS  THE  SUTHERLAND  VISCOSITY  LAW. 

RETURN 

1 RHO«1.224999*Y*(-.ll76033E-3*Y*< .43371 9E-04Y*(-. 746 1659E- 13 

1 ( .55376o3E-l8-.9572T27E-24*Y)T) J “ 

T«(1.0317O2E9-4.1O3O03E4*Y)/O 
■~G0  TO  0 

2 RHO-1.9901424Y* (-.2940114E-34Y*<.1993974E-74Y* (-.7637263E-12 
1 ♦Y*(.16i592iE-l6-.l476764E-2l»Y) ) ) ) 

T-216.65  __ 

GO  TO  8 

3 RH0«1.8156l4Y*(-.235749E-34Y*( .130807E-74Y*(-.3819651E-12 
“T^Y*  ( .§798720E-T7-.3656654E-02»Yn  )1 

T* ( 1*250 O50E9 46 .55341 6E3#Y) /D 
GD  TO  0 " 

4 RHO-1.  109444 Y*(-.l 140029E-34Y« (• 40 1740 1E-04Y«( -.1 039241 E-l 2 
1 ♦Y*dll38793E-17-.5052l35E-23*Y) ) D 

T« (0.039O03E04l.7938E4*Y)/O  

GO  TO  0 

5 RHO*.0974979E-l4Y« ( - .4 1 7905E-54 Y* ( . 3529753E- 1 0 ♦ Y* ( . 1 1 77 1 44E- 
1 4Y*(-.2567072E-194.l449ll3E-24*Y)  ) ) ) 

T“270 . 65 
GO  TO  0 

6 RHO«.1O29O02E-1^Y*(.1O01053E-54Y* (-.052361 9E-1 04 Y*( .2075003E 
1 ♦Y*(-.2l64824^i^.06O425E-?5lfYh  )1 

T*(2.301562E9-1.233088E4#Y)/D 
GO  TD  8 

7 RHO«0.4636»EXP(-0.12207E-3*Y) 

T«(3.157O00E9-2.4^3O41E4«Y)/O 
GO  TO  0 

END 

SUBROUTINE  FLIGHT < TIMEtU tKUTTA)  

DIMENSION  U<12) 

DIMENSION  TMACH(ll) .TKA(ll) t TKDYAW (11)  » TKL ( 1 1 ) tTKM(ll) . 

1 TKF(ll)  tTKTdl)  tTKHdll  tTKS(ll)  tTcPdll 

OAT A RE/6. 370E6/t OMEGA/O. 72915E-4/fPlOFOR/. 785390 l/»TWOG/l 9. 

~"re«nuwtn"at  Radius  df  tht  ea*th~at  thfto xtstor  in  heter^ 

OMEGA-ANGULAR  VELOCITY  OF  THE  EARTH  IN  RADIANS/SEC 

TABLE  OF  EQUIVALENCES 
U(l)  - X 
U ( 2 ) - Y 

\JT5T~li~  z " 


00046700 
00046000 " 
00046900 


00047000 
00047100 
00047200 
00047300 
~00O~47400“ 
00047500 
00 047600 
00047700 
000  47500 
00047900 
0 00400  00“ 
00040100 
00048200 
00048300 
00048400 
00040500 


14 


00048600 

00040700 


-14 


00040000 
00048900 
00 0490 00 
00049100 


OO049?OO 

00049300 

00049400 

00049500 

“00040600 

00049700 


"00040500- 

00049900 


00 05000 0 

00050100 

“00050200“ 

50418/00050300 

00050400“ 

00050500 
00050600 
00050700 
“00050800“ 
00050900 
00051000 


U ( 4 ) ■ SPIN 
U ( 5 ) ■ XDOT 
U ( 6 ) ■ YOOT 
U ( 7 ) ■ 200T 
U(0)  ■ SPIND 

C U ( 9 ) ■~XD8L 

C U ( 1 0 ) ■ YOBL 

C Udl)  ■ ZDBL 


00051100 

00051200 

00051300 

00057400 

00051500 

0005T500 

00051700 
00051500  “ 
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c 

IM 12)  « DUMMY 

00051900 

— 

EXTERNAL  CORAS  TJ  00520710 

COMMON  EMO»EMB*SPI*FC*BRATe*DELT*To»TB*ISW»V»THETA»FFCTR»CALSQ*  00052100 

lVWO.VXW.ALT.R*IEN0fCMACH.REYNL0.RE5IS.CAL*DL0N6.10PTY,YAW.AM0Mf  0O052200 

2 BMOM*PSI*WTAREA*ISEP*NABLE*I GLIDE  00052300 

CDmMDN/COFCOM/XCG.SMARG.EM.THID.PRNU  »ALTRTM.CNATRM,GLIOE.EPSTHE  00052*00 
1 *STAFAC*YAWNU*TKAtTKOYAW»TKL*TKM»TKF»TKT*TKH*TKS»TCP*TMACH*NARTBLOOO 52500 

COMMON/rfINCOM/RWC.XWC 
IF(ISW.EO.O)  GO  TO  10 

00052500 

00052700 

IF(TImE.GT.TO*TB*0ElT)  00  TO  9 

call  burn(time*xmass* thrust) 

0G0S2B07 

00052900 

_c 

eh.xmass- 

THRUST- I NOUCEO  DRAG 

00053000 

00053100 

FFC»FF7TR*TflT7~~  ' ' 

H*4. 44B23#THRUST 

53053200 

00053300 

TERmX«H«U15) 

TERMY«H*U(6) 

6o053*5ff 

00053500 

12  VSG«Ur5T<H»24UC6)**2VU(7)4*2 
V«SQRT(VSQ) 

UO0336W 

00053700 

WiKF+VTZT  ' 

XSQ«UU)**2 

00053800 

00053900 

YSG«UP##Z 
ZSQ*U ( 3 ) #*2 

00054000 

00054100 

c 

~R«SQRT rxSQ^YSQ^ZSQT 

”77054277”** 

00054300 

AUT-R-Rt  ~ 

00754400 

ALT  *U  ( 2 ) 4X50/2. /RE  0 0 054500 

0RC05X«U(1)/R  00054600 

ORCOSY«UP/R  00054700 

0RC0SZ«UT3T/R  ' ' * O0054B00  ~ 

CALL  SOUNO(AfGtRHOfVISCO)  « 00054900 


THIS  GENERATES  SPEEO  OF  SGUNO*GRAVlTY*AIR  OENSITY*  AND  VISCOSITY. 

I F fV.'EQ . 0 • 0 ) GO  TO  13 

termx«terhx/v 


NUMBERING  SKIN  FRICTION 


*07055777 

00055100 

0O055200 

00055300 


TERMT-TERMY/V 

••••COMPUTE  VELOCITY  OF  WIND  AS  A FUNCTION  OF  ALTITUOE. 

00755400 

00055500 

“C” 

HARG»1  .707*11  f27“ 

00055600 

C 

VW»VWO«RWC*VWC (HARG) 

00055700 

~C 

VCW»VXV4XVC*VWCiHARGr— “ 

v w*v wo 

00055BQ0  “ 

00055900 

vcw*vxw 

VRELSQ«  < <U<5> ♦VW)w»2^U(6)«*2^ (U(7) ♦VCW)«*2) 

00055700 

00056100 

VREL  *S  OR  T1  VREL5TJT”  4 

cmach-vrel/a 

0 0 05620  0 
00056300 

FRICT«0.0 

00056500 

IF(OLONG.EQ.O.O)  GO  TO  *8 

reynld»dldng*a#rhd#cmach/v jsco 

07056677 

00056700 

ALR-AL0G10TREYNCDT  ' 

00056000 

PWR«0.05*ALR 

00056900 

SFC*0.455/ALR#*2.5B/Ur40*2*CMACH*#27»*PWR~ 

00057000 
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IFUSEP.EO.O)  GO  TO  48 

F R I C T-WTAP  E A*  5F  C 

C DENOM-MASS  OF  PROJECTILE  IN  KG. 

48  DENDM-0 .4536#EM 

QYNPRS-0  *5#RHO#VRELSQ 
IF ( lOPTY.EQ. 1 ) GO  TO  20 
YAW-0.0 

YAWSQ-6Y0 

XKOYAW-O . 0 

xTerms-o.o 

yterms-o.o 

ZTERMS-0.0 

U ( 8 ) *0 . 0 _ 

2T  t all^CdraGI cm acH t ORaG t CaDT" 

CDFDRG«FFC*DRAG^XKDYAW#YAWsQ4FRICT4CAO 

this  GENERATES  THE  CORRECTED  COEFFICIENT  OF  DRAG 

DIFFERENTIAL 'EQUATIONS 

22  CONTINUE  _ 

form«piofdr*calsq*Oynprs 

DRG— COFDRG^FORM 
C 

C RESIS  IS  AIR  RESISTANCE  IN  POUNDS. 

RESIS-DRG  /A. 44823 

C***#  PROVISIONAL  CONSTANT  GLIDE  ANGLE  JAN  75 
rf( TGL! DE.NE.  1 ) GCTTD  150  " “ 


00057100 

00057200 

00057300 

00057400 

00057500 

~00O57633 

00057700 

60357803 

00057900 

OO0580OO 

00058100 

50058201} 

00058300 

8o358To3 

00058500 

00058600 

00058700 

00058800' 

00058900  _ 

00059000 

00059100 

00359200 

00059300 

00559400 

00059500 

55059650 


THE«ATAN(U(6)/U(5>) 

IF ( NA8LE «EQ • 1 ) GO  TO  62 
IF (THE. LE. GLIDE)  NA8LE-1 

herr-abs(the-glide; 

IF(HERR.LE.EPSTHE)  NA8LE-1 

IF(NABrE.NE.l)  3*TTo  60 

62  SAVE*0.5*ARSIN(DENOM*TWOG*COS (GLIDE)/FORM/CNATRM) 

ALPH-AM INI (ALTRlMf SAVE) 

SAVE-CNATRM*FORM 

U (4) »SaVE*SIN ( ALPH) /DtNOM 

c***»  U(4)  IS  COMPUTEO  AS  NORMAL  ACCELERATION  <M/S**2>  INSTEAD  OF  SPIN 
YAW-AlPlJ 

FL-0.5*SAVE#SIN(H.#ALPH) 

FDI*~SaVE*($IN(ALPH) >**2 
S INTH-U ( 6 ) /V 
CDSTH-U (5) /V 

termx«termx-fl*sinth*fdi*cdsth 
TE  R m Y - t E RM  Y ♦FT*  C35ThTF15T  5 TnTm  


00059700 
'06059863 
00059900 
"0 0063 060 
00060100 
00060200 
00060300 
06060405 
00060500 
06360600 
00060700 
83560800 
00060900 
000610OO' 
00061100 
06061200 
00061300 
60061435 


^60  CONTINUE  00061500 

PROVISIONAL  CONSTANT  GLIDE  ANGLE  JAN  75  ' 00061600' 

U(10)-(DRG»U(6)/VREL^TERMY)/DENOM-G#DRCOSY^0.53166E-8*UP  00061700 

1 ♦2.*OMEGA*U<7) ♦YTERMS  60  061833 

U ( 9 ) - ( DRG* ( U ( 5) ♦ VW ) / VREL  ♦TERMX)/DENOM-G*DRCOSX^XTERMS  00061900 

UTl  1 F--2.  *0HET^#UT6r*G*DRC65Z#ZTERM5VDRG*TlJ (T)  ♦VCWT/VREL/DENOM  5308'206T 
U ( 12 ) -0 • 0 00062100 

C AC  IS  THE  ACCELERATION  DF  THE  PROJECTILE  ALONG  THE  TRAJECTORY.  00062200 
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c 

AC*(U<5)*U<9>  *U(6)*U(10)  *U(7)*U(  11M/V 

00062300 

IFTIOPTY.'EQ.OT  60  fO  14 
U(8)»-RH0*CALSQ**2/AM0M*XKA*U<4) *V 

000624OO 

00062500 

14  IF<KUTTA.EQ.4J  ThETa»AR5IN <U (6) /V ) *57. 3 
RETURN 

06062600 

00062700 

13  U<9)»0. 

U < 1 0 ) *-G 

0O062BO0 

00062900 

THTETA*Wi 

RETURN 

00063000 

00063100 

9 EM*EM0 
GO  TO  11 

000«2H50 

00063300 

10  EM*EMO 

11  TERMX*0 • 

Oo 063460  " 
00063500 

TERMY40. 

FFC«FFCTR 

00063600 

00063700 

60  TO  12 

20  CALL  AC0EFS<CMACH.YA**XKA*XKDYAW.XKL.XKM.XKF.XKT. 

00063&00 

00063900 

lXKH.XhS.XCPt 
VXRL“U (5) *VW 

00064000 

00064100 

— v7RL*u(7y*vctr 

VRLSO* VXRL**2*U (6) **2*VZRL**2 

00064200 

00064300 

i 

c 

VRL*S0RT (VRlSO) 

TEST  FOR  DYNAMIC  STABILITY. 

UO064400 

00064500 

BOTTOM»87D#BMOH*DYNPRS*CAl.*43*XKM" 
STAFAC*  <U (4) *AMOM) **2/B0TT0M 

D0T)6460D 

00064700 

Tr  (5T A pACTLeTITo)  001015 
C****  COMPUTE  THE  YAWING  FREQUENCY 

TTOO640OO 

00064900 

YAWNU*AM0M/BM0M»SQRTU.-r./BTAFACr#U(Al/6.2632 
C###«  COMPUTE  THE  OVERTURNING  MOMENT  AND  PRECESSIONAL  FREQUENCY 

W065OOi) 

00065100 

0TNM0M»2  , ♦XKW«CXL"3#DYNPKS 
PRNU*OTNMOM/AMOM/U (4) /6.2032 

WO 652 W 
00065300 

C 

“COMPUTE-  TAW  OT~REPOSE.  - ' 

ALPHAB«RH0*CAL*VRLSQ*<XKL*XKM*VRLSQ*CALSQ*XKF*XKT*U<4>**2> 

“00065400 

00065500 

IFTABS  ( ALPHABt  VLT. I.E-20)  GO  TU‘25  ~ *“ 

ALPHAA«AM0M*XKL*U(4)/CALSQ/ALPHAB 

D0065600 

00065700 

“ACPHAB«DEN0W#XKT#Uf41  /ALFHAB 
AMB-ALPHAB-ALPHAA 

U0065800 

00065900 

ALFHAX*AHB#  (UT6T«UC1  IT-YTRL^UTTO  J 1 - ALPHA  B*VZRl*TT 
ALPHAY* AMB*  ( VZRL*U  ( 9 ) -VXRL*U  (ID) 

00066000 

00066100 

ALPHAZ»AMB*  ( VXRL*U TI 0)  -UC6 J *U  191  ) ♦ ALPHAB*VXR1_#G 
YAWS0*ALPHAX**2* ALPHA Y**2*ALPHAZ**2 

0OO66ZW 

00066300 

YAW«SORT(YAWSQ) 

IFIYAW.GT. 1.5708)  60  TO  25 

WO  66455 
00066500 

AlTB*  (VXRL*ALPHAZ»VZRL#ACPHaXI  *vrl 

AR61*(VXRL*ALPHAY-U(6) *ALPHAX)*VXRl-(U(6)*ALPHAZ-VZRL*ALPHAY) 

00066600 

•VZRL00066700 

J ' IF (ABS(AROl) .LE.1.0E-20)  Go  TO  50  - 

PSI»57.3*ATAN(ARG/ARG1) 

50B66WO 

00066900 

GO  TO  53 

50  IF ( ARG#ARG1 ) 51*52t52 

TT0O67OOTJ 

00067100 

5I~P5T*  -9  0 V 
GO  TO  53 

00067200 

00067300 

"52  PST*90  • 

T)00674W~ 

-109- 


o on  oj  o 1 olooooo 


53  CONTINUE 
C 

PSI» ORIENTATION  OF  YAW.  THIS  IS  THE  ANGLE  BETWEEN  THE  PLANE 
CONTAINING  BOTH  THE  VELOCITY  aKIO  Ya'N  VECTORS  AND  A VERTICAL 
PLANE  CONTAINING  THE  VELOCITY  VECTOR.  IT  IS  MEASUREO 

clockwise  from  the  vertical  plane. 

END' OF  COMPUTATION  oftaw. ; 

OKFN»CAL*XKF*U (A)  

xklvSq»XkL*vrlso 

RDSO*RHO*C ALSO/OENOM 

XTERMS-RDSQ* (XkLVSO*ALPHAX*DKFN* ( AlPHAY*VZRl-ALPHAZ*U <6) T7 
YTERMS*R0S0*(XKLVS0*ALPHAY*DKFN*(ALPHAZ*VXRL-ALPHAX*VZRL) ) 
ZTERMS.P05O*  (XKLVEO»ILPNJZIDKrN*  ( AlPHAX*uT6 f-^CPHAY *VXRtTT 

XKDYAW«2.54647*XKDYAW  _ 

IF (YAW, LT. 0.69)  GOTO  21 
~ PRINT  55»RffO*XTERMS,YTERM5,ZTERMT 
55  FORMAT ( 1H  *1P4E10.5) 

•*••*«*••*•*••*•• 4 ••*•••* 

GO  TO  21 

'25  PRINT  26.STAFAC 

26  FORMAT <1H0»40HUNSTA8LE  PROJECTILE  STABILITY  FACTOR  ■ .F10.4) 
CALL  EXIT' 

END 

SUBROUTINE  CCRaS fEffACN.DRAS'tCAin 

program  computes'tne  coefficient  of  drag  versus  mach  number 

and  the  THE  COEFFICIENT  INCREMENT  due  to  canaros. 


00067500 

0O067666 

00067700 
00067800 
00067900 
06068660 
00068100 
' 66668260 
00068300 
'00068400 
00068500 
00068600 
00068700 

OO068S0S 

00068900 
' 0 0 06  90  0 0 
00069100 
00069200 

00069300 

00069400 
00069500 
0 00696 00 
00069700 
00066806 
00069900 

S0O7SSSS 

00070100 

00070200 

00070300 

60670666 


OIMENSION  XMTBL(ll) *C0T8L<11) *COT8l(11)  00070500 

COMMON  EMO.EMB , SPI  ,7Ct  BP  ATE, CELT , T DTTH , T5¥7v  .THETa.FFCTR.CaLSQ, 0o0?0466 

1 VW,VCW.ALT.R.IENO«CMACH,ReYNLO*RESIS*CAL*OLONGi IOPTY.YAW.AMOM,  00070700 

2 BMOM.PSI * WTaREA. ISEPVNaBLe* iGLfDE  00076805 

COMMON/ORGCOM/  XMT8L * CDT8L . COT8L *NT8L  00070900 

IF(NTBl.NE.O)  GO  TO  5 ""  " " O0OTI6OO 


IF(EMACH.LE.0.80)  go  to  1 

TFIEMACH.LE.I.IDT  50'TITl 

IFIEMACH.LE.3.0)  GO  TO  4 
EM3«EMaCH-3.0 

DRAG»0.09*EM3*(“0.02*0.002*EM3) 

RETURN 

1 ORAG»0.0589 
RETURN- 


00071100 

06071200 

00071300 

066TH66 

00071500 

'600716*00 

00071700 

06071860 


3 C*10«*(EMACH"0.8)  00071900 

OR AG«0. 07736 *C**3*EXP(“C)*0. 0589  00672666 

RETURN  00072100 

4 DR AG«0 • 21547*E MACH* ( “0,05134* 0.003 i7*EMACH)  06072208 

5 00  6 J-l.NTBL  00072300 

ITTTmaCN . IT .THTBUrJi  ySO  T6“8 66072A60 

6 CONTINUE  00072500 

8 JL»J»1  " ~ " 00672600 
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NFV  HASTEN 


~TEBUPDTE“tOG  PABFWIS 


FRAC- (EMACH-XMTBL ( JL) ) / ( XMTBL  < J) -XMTBL ( JL) ) 

CD-CUTBL7JIT)  ♦TCDTBL  (J)  -COTBUrOLT)  *frac 

CAO-O.O 

~TFTN  ABlE^F.TTWTBTO  “ 

CAO-COTBL (JL) ♦ (COTBL ( J) -COTBL  < JL) ) *FR  AC 

TO  CONTINUE  ~ 1 

ORAG-CO 
CAD-CAO 


00072700 


TOO  728  00 

00072900 
“00073000  ~~ 
00073100 
"Oo  07  3200 
00073300 

~oirm*oo 

00073500 


RETURN 

ENO  0oOT35OO 

SUBROUTINE  ACOEFS (EMACH# Y AW* XKA*XKOYAW*XKLtXKH*XKFtXKTtXKH*XKSfXCPOO 073700 

IT 


000738W 

00073900 

00074000 


DIMENSION  TMACH(ll) tTKA(ll)  » TKDYAW (11)  • TKL (1 1 ) t TKM ( 1 1 ) t TKF ( 1 1 ) , 

T TK T (11 ) *TKH  < 11  )'tTRS  ( 1TFYT CFTTll 
COMMON/COFCOM/XCG* SMARG tEM t THI 0*  PRNU  # ALTRI M# CNATRM* GLI OE tEPSTHE  00074100 
T * S T A F A C7  Y A W NU  t TK  A*  TKD  YA  WVT  KL  VT  KM  t TKF  FT  KT  7TK  H rTKS  * T C P t T MAC  H t~NA  PTGUo  00  74T0OT 
EMACH  • MACH  NUMBER  00074300 

XKa'-'SPIN'DaUPING 'MOMENT  COEFFICIENT  00074*00 

XKOYAW  ■ YAW  DRAG  COEFFICIENT  00074500 


TKL'-  LIFT  FORCE  COEFFICIENT" 

XKM  « OVERTURNING  MOMENT  COEFFICIENT 
Tier-  MAGNUS  TOBCF  COEFFICIENT 
XKT  ■ MAGNUS  MOMENT  COEFFICIENT 
TKRT-  DAMPING  MOMENT  COEFFICIENT 
XKS  - PITCHING  FORCE  COEFFICIENT 


“06T7T5W" 
00074700 
T>  0074800' 
00074900 
“00075000 
00075100 


XCP  « CENTER  OF  PRESSURE  AFT  OF  NOSE  ^“CALIBERS 


C 

c 

C~ 

c 

<T“ 


RELATIVE  TO  T387  TYPE  PROJECTILE. 
XKT-0.0 


MEMO.  RPT.  NO.  2023“ 


XKFS0.  157 

SYAW»SIN(YAW)**2 

GOTO  50  — 

51  CONTINUE 

~ ^BOTTTJ-ltNARTBL  ~ 

IF<EMACH.LT.TMACH(J) ) GO  TO  70 


00075200 
00075300 
~00075*00 
00075500 
U 0075500 
00075700 
*00075500 
00075900 
00076000 
00076100 


“60  CONTINUE  

70  JL-J-1 

FfUC-TENACH-TMACHrJLTT/(TMACHTJT-TMACM(JL7T 

IF(TKA<J) .EO.O.O)  Go  To  52 

XKAT-  TKATJLT ♦ fTKA  TUI  -TKA  TULTPFWAC 

52  IF  (TKOYAW(J) .EQ.0.0)  GO  TO  53 


XKDTAW»YKDTA¥(  JL)  ♦TkDTAWI  Jj-TKDTAWl  JLl)  R AC 

53  IF  (TKL(J) .EQ.0.0)  GO  TO  54 

XXL— TKL  f JL7  ♦ TTKL  ( JT •TKLTJL)  7 *FR  A C 

54  IF(TKM(J) .EQ.0.0)  GO  TO  55 

XKR*TKH  TJLT^  (TKHTJ)  -TKM  ( JL)  7 -FRAC 

55  IF(TKF(J) .EQ.0.0)  GO  TO  56 


005752075 
00076300 
00076400 
00076500 
“00076600 
00076700 
*00076800 
00076900 
00077000 
00077100 
00077200 
00077300 
"000  7 7 AW 
00077500 
00077600 
00077700 
”00077800 


XKF-TKF  ( JLTV  {TKF (J7- TKF  C3L')  I-FKAC 

56  IF(TKT<J) .EQ.0.0)  GO  TO  57 
XKT-TKTTJL)  ♦TTKTIJ)  -I 
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57 

IF(TKH(J) .EQ.O.O)  00  TO  58 

00077900 

xkH»tKH(JL)* (tkhcji-tkhiJL) )*Frac 

00078000 

58 

IFITKSIJ) .EQ.O.O)  60  TO  59 

00078100 

XKS-TKS(JL) ♦ (Tk5( J7-TK51JL) )*fRaC 

00078200 

59 

IF ( TKM ( J) *EQ • 0 • 0 ) Go  TO  62 

00078300 

iFixkl. eq.o.o)  Call  exit 

00078400 

SMARG— XKM/XKL 

00078500 

XCPmXCg*S"MARG 

ooOtAGbo 

RETURN 

00070700 

62 

IFITCPIJ) .EQ.O.O)  SO  TO  63 

00078800 

XCP-TCP(JL) ♦ (TCP ( J) “TCP ( JL) )*FRAC 

00078900 

63 

SMARS-XCP-XCG 

00079000 

xkm«-xkl*smarg 

00079100 

return 

60079200 

50 

CONTINUE 

00079300 

IF (EMACH.LE.0.8)  GO  TO  10 

00079400 

IFIEMACH.LE.0.9)  GO  TO  20 

00079500 

IF(EHACH.LE.I.O)  GO  TO  30 

OO0T96OO 

IF(EHACH.LE.l.l)  go  TO  35 

00079700 

1E(eha'ch.le.i.30)~^oTO  4® 

00079800 

IF (EMACH.GT. 1 .5)  GO  TO  45 

00079900 

c 

VALID  for  emaCh  GTR  0.8~aNd  Lt  1.5 

00080000 

5 

XKA»0.0038*0.002*EXP(-1.5*<EMACH-0.8) ) 

00080100 

c 

VALID  FOR  ENACT  GTR  0.0 

00080200 

6 

EM9-EMACH-0.9 

00080300 

~H  OLD  - r.  -EXP  f -5.»EM9  y 

0OO8O4OO 

XKL»0. 5507*0. 4*H0L0 

00080500 

XKL»XkL*6.6*SYa! 

OOO8O600 

- 

XCP-O. 237*1. 57*H0LD 

00080700 

c 

VALID  rOR  EMaCm  gtr  itd 

00000000 

7 

XKS«-4#0*1 *78*<EMACH-1 •) 

00000900 

c 

ViLTD  FOR  EHACTTSTR  T7T 

00061000 

XKDYAW-1 .5*2.38*EXP (-2.72* (EMACH-l .1 ) ) 

00081100 

“ c 

valid  for  emach  gtr~it3 

00081200 

XKH-3.7 

00081300 

"C 

VALID  FOR  ALL  E^aCH 

00001*00 

9 

smarg»xcp-xcg 

00081500 

XiOi— TCkL^MaRS  ■ ■ 

00081600 

IF(NARTBL.NE.O)  GO  TO  51 

00081700 

RETURN' 

OOD&IBOO 

10 

XKA-0.0058 

00081900 

XKDYAW»1 *5 

60Of2O00 

11 

XKLW0*62“0* 077 #E MACH 

00082100 

T(KLiX)<L*4  .3*5Ya¥ 

0008?2O® 

XCP»1.2-1.07*EMACH 

00082300 

12 

XKS--4.0 

00002400 

13 

XKH*0*71*2*3#EMACH 

00062500 

60  TO  $ 

00082600 

20 

EM8*EMACH-0.8 

00082700 

"XKiVT r.  orr3B*rr,oo2*Exp(-i  .^ehri 

00082801) 

XK0YAW>1.5«2,5*SIN(6.283*Em8) 

00082900 

730  TO  T1  ' 00  0*530  Off 
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"NEtT  MASTER 


'7EBUPDTE~U)G  FAGE"  ^017 


30 

EM8*EMACH-0 • 8 

00083100 

XK A*0.0038*0.002*EXP  l " 1 . bwEM 0 T 
XKDYAW»1.5*2.5*SIN(6.283*Em8> 

00003200 

00083300 

EM9»EMaCB-0,9 

H0L0»1.-EXP(-5.*EM9) 

00083400 

00083500 

XKL-0. 5507*0. 4»HOLO 
XK|_«XKL*5.5*SYA» 

00083600 

00083700 

■XCP.0 . ZT37  ♦T','57»H0LG 
GO  TO  12 

"00003800 

00083900 

35 

EM8*EMACH-0 • 8 

XKA*0.0038^0.002*EXP(-1#5#E:M8) 

00084000 

00084100 

XK0YAW»1.5*2.5*SINt6.283*EM8) 

EM9-EMACH-0.9 

00084200 

00084300 

>TOEtr-l . -EXP  t *EH9) 

XKL*0.5507*0.4*HDLO 

TJUOffSSTTO 

00084500 

XKL«XKL*5.$*$rAW 

XCP«0*237M«57*HOLD 

00084605 

00084700 

XKS*-5.78^1  •78*EMACh 
GO  TO  13 

00084000 

00084900 

40 

EM4EMAEH-0.8 

XKA*0.0038^0.002*EXP(-1 .5#EM8) 

00085000 

00085100 

XKDYa*»»I. 5*2. 38*£XPT-2. 72*  (EMaCH- l.TD 
EM9»EHACH-0.9 

0008520c 

00085300 

HDLD*1 •-EXPf-b#*EH9) 
XKL*0*5507*0«  4#H0L0 

00085400 

00085500 

XKL«xkl*6.6*SYAW 
XCP-0. 237*1. 57*HDL0 

00085600 

00085700 

XKS— 5.78*1 ,78*EMACtf 
GO  70  13 

00085800" 

00085900 

4b  XKA*0»  003B*0  • 002*EXP (-T75w  (EMACH-0  .U)  1 
GO  TO  6 

00086000 

00086100 

9 

END 

SUBROUTINE  RUNGE1 (V.WtNEQtNDRDtDIFEQ) 

~OOOB6200 

00086300 

DIMENSION  V(I2)fWI48j 
NV«NCQ*NORD 

~W08540U 

00086500 

N*NV*NEO 

RETURN 

0008660TT 

00086700 

“C 

ENTRY  RUNGE2 ( T t DT ) 

~0OO860QO 

00086900 

DT2«0rr*,5 

DT6«DT/6. 

GOOB70GG 

00087100 

1 

otniir(N  

W(I)«V(I) 

oo0872u0 

00087300 

DO  2 JSTf3~ 
NJM*NEQ^N*(J-1) 

OQUB  7*UU 
00087500 

JDECK»N*J 
IF  ( J”3) 3***4 

00087600 

00087700 

3T 

DTW-DT2 
GO  TO  5 

00087800 

00087900 

> 5 

DTW-DT 

TW«T»OTW 

OO  0060  0 0 
00008100 

vxr  6 r*iiNv 

00008200 
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orvoninooooniooo 


K* I ♦ JDECK  , 

— L-TN3W : 

W<K)*W<I)*W<L>«DTW 

6VTT)«W(K) 

call  DlFEQ(TW#V#J) 

~~D0  2 1*1 tH  ” 

K* I ♦ JOECK 

2rvrKr*vm  * 

00  7 I* 1 #NV 
Kl*MNE0 
K2*KUN 
K3*K2*N“ 

K4*K3*N 

rm  )*lT(Tr^T^W7m^2^'(VTX21  iTrnr3lTV»T>01T 

T*TW 

CALL  DlFEQ  <T # V#4)  ' 

RETURN 
“ "END 


00066300 

0OO984OO 

00088500 

00088607$ 

00088700 

66668366 

00088900 

00089000 

00089100 

00389200 

00089300 

00089400 

00089500 

60089663 

00089700 

66689800 

00089900 

00090006 


FUNCTION  VWC(Y) 

IF ( Y*OE*6700* ) 60T03 

VWC-5.181644,972E-5«YM,3494E-7*Y«Y 

RETURN 

3 VWC*l0.058O.962MC0S  ( *42946E-3*  (Y-9448.8)  ) 

IF ( Y»GT • 13700* ) WRlTEt6#IT 
RETURN 

1 FOR H aT7 28TT ALTITUDE"  ABOVE  ~l37WHrr£A5$ 

ENO 

SUBROUTINE  BURN (TIME# XNASSfTHRUSTJ 

SUBROUTINE  COMPUTES  PROJECTILE  MASS  IN  POUNOS  MASS  ANO 
ROCKET  THRUST  IN  FOUNDS  FORCE  ~ 

TO-TIME  AT  WHICH  BURNING  COMMENCES 

EMO*  I N HTAL"  74A5S  #~LBM 

EMB«BURNT  MASS # L8M 
TIME*TIME  AFTER  LAUNCH  #SEc 
SP I*SPEC IF IC  IMPULSE#  L»F/LBM/SEC 
“ FC-CONSTANT  NOMINAL  THRUST  LEVEL#  LBF 


00090100 

00090200 

00090300 

00090400 

00090500 

66690660 

00090700 

66696S66 

00090900 
63  6917)66 
00091100 
66091260 
00091300 
00091400 
00091500 
067591660 
00091700 
66691660 


IBURN*  INDICATOR  OF  COMMENCEMENT  OF  BURNING ( I BURN* 1 ) 

’ UECT-F15F  TIME  OF  THRUST^  ASSUMED" 


00091900 


BRATE*FC/SPl "BURNING  RATE#  LBM/SEC  00092100 

TB*(EMO-EMB)/BRATE*EFFECTTvE  BURNING  TIME#  SEC  66692260 


00092300 

COMMON  EM0»EMB#$PI#FC#8RATE#DELT  #TO*TB* ISW# v#THETa»FFCTR#CaLSQ»  “0757)92460 


1 VW,VCW#ALT#R*IEND#CMACH#REYNL0#RESIS»CAL*DL0NG#I0PTY#YAW#AM0M,  00092500 

2 1JM6H7P51  # wTaRE* #T5EF7NaBLE # T3LTDE~  “ 1 66692600 

IF(TIME.LE.TO)  GO  TO  1 _00092700 

IF(TIME.LE.TO^OELT)  GO  TO  2 06692666 

T2*TO^TB  00092900 

IF(TIME.LE.T25  GO  TO  6 “6609T066 

IF(TIME.LT,T2^0ELT)  GO  TO  4 00093100 

’ THRUST  4 or; 06093260 


XMASS*EMB  00093300 

RETURN  00693466 


* 


« 
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NEW  MASTER 


IEBUPDTE  LUG  PAGE  U019 


1 xmass«emo 

mRUSTmO. 

RETURN 

2 XMASS»EMO-BRATE# (TIME*T0)#*2/(H«#DELTJ 
THRUST* (TIME* TO) /DELT^FC 

RETURN 

3 xma$s*emo-brate*<time-to-delt/2.> 
THRUST.Fd  ~ 

RETURN 

4 XmaSS«EMO-BRaTE* < < TB-DELT/2. ) ♦ <TTME-T2 > * { 1 
THRUST*FC*(1.-(TIHE.T2)/DELT) 

RETURN 
END  * 

FUNCTTON  IBURNmMErALT#QEfVELOT 


00093500 

TTDTX93  BUG— 

00093700 
00093600 
____  00093900 

00094000 
00094100 
00M4T015 
00094300 

(TTME-T2) /DElT/2.  F)  OO0944O0 

00094500 

U0D94600 

00094700 

0OT9'VBU0 

00094900 


FUNCTION  PRODUCES  INDICATION  OF  COMMENCEMENT  OF 
BURNING--IBURN«U 

C — I BURN*0~~ UNT IL  BURN  BEGINS 

C USER  CHOOSES  TO  FROM  CURRENT  TIME  OR  ALT  i ALTITUDE)  OR 

C QFTLUCAC  OGATTRANTTLETaTTOTTP  OR  VELtPTVELOCITY) 

COMMON  /SRCOM/TM 

DATA  ALTMAX/30000.7%WMIN/0,0/~ 

IF(TIME.GE,TM)  GO  TO  3 
IFULT.GE.ALTMAX)  ~G 0 TO  ^ 

IF(VELO.LE.VMIN)  GO  TO  3 

u rrfor;DT,45,u)  go"to~2 ~ 

IBURN-0 

RETURN 

C 2 IF ( QE • LT • 46 . 0 ) GO  TO  3 
C 4 ~IBURN»0 
C RETURN 

P I BURR*  I 

RETURN 

END 


00095000 

00095100 

06095200 

00095300 

0oS9S4U&~" 

00095500 
00095690 
00095700 
00095000 
00095900 
OTT09600T5 
00096100 
00096200 
00096300 
00096400 
00096500 
U0GR66U0 
00096700 
00096B0  0 


TN  NM  07  RECTORY.  TTR  TS  NOW  ALTERED . 
0000000 
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AERU  Data  (HRL  CALC)  for  M509  8 IN  lew  PROJECTILE 


M*CH  no 

COEF  DRAG 

DRAG  I N C R 

0*0 

0.1300 

0.0 

O.T500 

■OVT30D 

C.IT 

0.8500 

0.1400 

0.0 

0.9000 

6.1550 

0.0 

1 .0000 

0.3000 

0.0 

1.05015 

0 .3800 

O.tf 

1.1000 

- 0.3600 

0.0 

1 ; 50  63 0.^3  iTo  dvo 

2.0000 0*2740  0*0 

2.5000  0.2390  " 0.0 

EWaBlE  TTMe  * 100.0000  THRUST  DRAD  FACTOR  * 


1.0000 


air  dens  factor  * ^i;oooo 


* 
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